$IBFTC TANDEM  DEBUG
      COMMON SRW,ITER,IEND,LER(2),NER(2)
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /GEOMIN/ CHORD(4),STGR(4),MLE(4),THLE(4),RMI(4),RMO(4),
     1RI(4),RO(4),BETI(4),BETO(4),NSPI(4),MSP(50,4),THSP(50,4)
      COMMON /RHOS/ RHOHB(100,4),RHOVB(100,4)
      COMMON /BLCDCM/ EM(50,4),INIT(4)
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      CALL TIME1(T1)
   10 IEND = -1
      ITER = 0
      DO 20 SURF=1,4
   20 INIT(SURF) = 0
      CALL INPUT
      CALL PRECAL
   30 CALL COEF
      CALL SOR
      CALL TIME1(T2)
      TIME= (T2-T1)/3600.
      WRITE(6,1000) TIME
      CALL SLAX
      CALL TANG
      CALL VELOCY
      CALL TIME1(T2)
      TIME= (T2-T1)/3600.
      WRITE(6,1000) TIME
      IF(NER(2).GT.0) GO TO 10
      IF (IEND) 30,30,10
 1000 FORMAT (8HLTIME = ,F7.4,5H MIN.)
      END
$IBFTC INPUT   DEBUG
      SUBROUTINE INPUT
C
C  INPUT READS AND PRINTS ALL INPUT DATA CARDS AND CALCULATES HORIZONTAL
C  SPACING (MV ARRAY)
C
      COMMON SRW,ITER,IEND,LER(2),NER(2)
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /GEOMIN/ CHORD(4),STGR(4),MLE(4),THLE(4),RMI(4),RMO(4),
     1RI(4),RO(4),BETI(4),BETO(4),NSPI(4),MSP(50,4),THSP(50,4)
      COMMON /RHOS/ RHOHB(100,4),RHOVB(100,4)
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
C
C  READ AND PRINT ALL INPUT DATA
C
      WRITE(6,1000)
      READ (5,1100)
      WRITE(6,1100)
      WRITE(6,1110)
      READ (5,1030) GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF
      WRITE(6,1040) GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF
      WRITE(6,1120)
      READ (5,1030)BETAI,BETAO,CHORD(1),STGR(1),CHORD(3),STGR(3),
     1MLE(3),THLE(3)
      WRITE(6,1040)BETAI,BETAO,CHORD(1),STGR(1),CHORD(3),STGR(3),
     1MLE(3),THLE(3)
      WRITE(6,1130)
      READ (5,1010) MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP
      WRITE(6,1010) MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP
      DO 10 J=1,4
      IF (J.EQ.1) WRITE(6,1140)
      IF (J.EQ.2) WRITE(6,1150)
      IF (J.EQ.3) WRITE(6,1160)
      IF (J.EQ.4) WRITE(6,1170)
      WRITE(6,1180) J,J,J,J,J
      READ (5,1030) RI(J),RO(J),BETI(J),BETO(J),SPLNO
      WRITE(6,1040) RI(J),RO(J),BETI(J),BETO(J),SPLNO
      NSPI(J)= SPLNO
      NSP = NSPI(J)
      WRITE(6,1190) J
      READ (5,1030) (MSP(I,J),I=1,NSP)
      WRITE(6,1040) (MSP(I,J),I=1,NSP)
      WRITE(6,1200) J
      READ (5,1030) (THSP(I,J),I=1,NSP)
   10 WRITE(6,1040) (THSP(I,J),I=1,NSP)
      WRITE(6,1210)
      READ (5,1030) (MR(I),I=1,NRSP)
      WRITE(6,1040) (MR(I),I=1,NRSP)
      WRITE(6,1220)
      READ (5,1030) (RMSP(I),I=1,NRSP)
      WRITE(6,1040) (RMSP(I),I=1,NRSP)
      WRITE(6,1230)
      READ (5,1030) (BESP(I),I=1,NRSP)
      WRITE(6,1040) (BESP(I),I=1,NRSP)
      WRITE(6,1240)
      READ (5,1010) BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      WRITE(6,1020) BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      IF (MM.LE.100.AND.NBBI.LE.50.AND.NRSP.LE.50.AND.NSPI(1).LE.50
     1.AND.NSPI(2).LE.50.AND.NSPI(3).LE.50.AND.NSPI(4).LE.50) GO TO 20
      WRITE (6,1250)
      STOP
C
C  CALCULATE MV ARRAY
C
   20 HM1 = CHORD(1)/FLOAT(MBO-MBI)
      IF(MBO.GT.MBI2.AND.MBI.NE.MBI2) HM1 = MLE(3)/FLOAT(MBI2-MBI)
      HM2 = 1.E30
      IF(MBI2.NE.MBO) HM2 = (MLE(3)-CHORD(1))/FLOAT(MBI2-MBO)
      HM3 = CHORD(3)/FLOAT(MBO2-MBI2)
      IF(MBO.GT.MBI2.AND.MBO.NE.MBO2) HM3=(CHORD(3)+MLE(3)-CHORD(1))/
     1FLOAT(MBO2-MBO)
      MBOT = MIN0(MBO,MBI2)
      CDMBOT = AMIN1(CHORD(1),MLE(3))
      DO 30 IM=1,MBOT
   30 MV(IM) = FLOAT(IM-MBI)*HM1
      MBIT = MAX0(MBO,MBI2)
      CDMBIT = AMAX1(CHORD(1),MLE(3))
      DO 40 IM=MBOT,MBIT
   40 MV(IM) = CDMBOT+FLOAT(IM-MBOT)*HM2
      DO 50 IM = MBIT,MM
   50 MV(IM) = CDMBIT+FLOAT(IM-MBIT)*HM3
C
C  CALCULATE MISCELLANEOUS CONSTANTS
C
      NER(1)=0
      NER(2)=0
      PITCH = 2.*3.1415927/FLOAT(NBL)
      HT= PITCH/FLOAT(NBBI)
      DTLR= HT/1000.
      DMLR = AMIN1(HM1,HM2,HM3)/1000.
      BV(1) = 0.
      BV(2) = 1.
      BV(3) = -WTFLSP/WTFL
      BV(4) = 1.+BV(3)
      MBIM1= MBI-1
      MBIP1= MBI+1
      MBOM1= MBO-1
      MBOP1= MBO+1
      MBI2M1= MBI2-1
      MBI2P1= MBI2+1
      MBO2M1= MBO2-1
      MBO2P1= MBO2+1
      MMM1 = MM-1
      CP = AR/(GAM-1.)*GAM
      EXPON= 1./(GAM-1.)
      TWW= 2.*OMEGA/WTFL
      CPTIP= 2.*CP*TIP
      TGROG= 2.*GAM*AR/(GAM+1.)
      CALL SPLINT(MR,RMSP,NRSP,MV,MM,RM,SAL)
      CALL SPLINT(MR,BESP,NRSP,MV,MM,BE,DBDM)
C
C  CALCULATE GEOMETRICAL CONSTANTS
C
      CHORD(2) = CHORD(1)
      CHORD(4) = CHORD(3)
      STGR(2) = STGR(1)
      STGR(4) = STGR(3)
      MLE(1) = 0.
      MLE(2) = 0.
      MLE(4) = MLE(3)
      THLE(1) = 0.
      THLE(2) = PITCH
      THLE(4) = PITCH+THLE(3)
      RMI(1) = RM(MBI)
      RMI(2) = RM(MBI)
      RMI(3) = RM(MBI2)
      RMI(4) = RM(MBI2)
      RMO(1) = RM(MBO)
      RMO(2) = RM(MBO)
      RMO(3) = RM(MBO2)
      RMO(4) = RM(MBO2)
C
C  INITIALIZE ARRAYS
C
      DO 60 I=1,2000
      U(I) = 1.
      K(I) = 0.
   60 RHO(I) = RHOIP
      DO 70 IM=1,100
      DO 70 SURF=1,4
      RHOHB(IM,SURF) = RHOIP
      RHOVB(IM,SURF) = RHOIP
   70 ITV(IM,SURF) = -10000
      RETURN
 1000 FORMAT (1H1)
 1010 FORMAT (16I5)
 1020 FORMAT (1X,16I7)
 1030 FORMAT (8F10.5)
 1040 FORMAT (1X,8G16.7)
 1100 FORMAT (80H
     1                         )
 1110 FORMAT (7X,3HGAM,14X,2HAR,13X,3HTIP,12X,5HRHOIP,12X,4HWTFL,11X,6HW
     1TFLSP,10X,5HOMEGA,12X,3HORF)
 1120 FORMAT (6X,5HBETAI,10X,5HBETAO,11X,6HCHORDF,11X,5HSTGRF,10X,
     16HCHORDR,10X,5HSTGRR,12X,4HMLER,11X,5HTHLER)
 1130 FORMAT (41H   MBI  MBO MBI2 MBO2  MM  NBBI  NBL NRSP)
 1140 FORMAT (53HL    BLADE SURFACE 1  --  UPPER SURFACE - FRONT BLADE)
 1150 FORMAT (53HL    BLADE SURFACE 2  --  LOWER SURFACE - FRONT BLADE)
 1160 FORMAT (52HL    BLADE SURFACE 3  --  UPPER SURFACE - REAR BLADE)
 1170 FORMAT (52HL    BLADE SURFACE 4  --  LOWER SURFACE - REAR BLADE)
 1180 FORMAT (7X,2HRI,I1,12X,2HRO,I1,12X,4HBETI,I1,11X,4HBETO,I1,11X,5HS
     1PLNO,I1)
 1190 FORMAT (7X,3HMSP,I1,2X,5HARRAY)
 1200 FORMAT (7X,4HTHSP,I1,2X,5HARRAY)
 1210 FORMAT (16HL      MR  ARRAY)
 1220 FORMAT (7X,11HRMSP  ARRAY)
 1230 FORMAT (7X,11HBESP  ARRAY)
 1240 FORMAT (52HL    BLDAT  AANDK  ERSOR  STRFN  SLCRD  INTVL  SURVL)
 1250 FORMAT (41H1 MM,NBBI,NRSP,OR SOME SPLNO IS TOO LARGE)
      END
$IBFTC PRECAL  DEBUG
      SUBROUTINE PRECAL
C
C  PRECAL CALCULATES ALL REQUIRED FIXED CONSTANTS
C
      COMMON SRW,ITER,IEND,LER(2),NER(2)
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      EXTERNAL BL1,BL2,BL3,BL4
C
C  CALCULATE LAMBDA AND VI
C
      BETAI = BETAI/57.295779
      BETAO = BETAO/57.295779
      TBI = SIN(BETAI)/COS(BETAI)
      TBO = SIN(BETAO)/COS(BETAO)
   10 RHOT = RHOIP
      RHOVI = WTFL/BE(MBI)/PITCH/COS(BETAI)/RM(MBI)
   20 VI = RHOVI/RHOT
      LAMBDA = RM(MBI)*(VI*SIN(BETAI)+OMEGA*RM(MBI))
      TTIP = 1.-(VI**2+2.*OMEGA*LAMBDA-(OMEGA*RM(MBI))**2)/CPTIP
      IF(TTIP.LE.0.) GO TO 30
      RHOMBI = RHOIP*TTIP**EXPON
      IF(ABS(RHOMBI-RHOT)/RHOIP.LT..000001) GO TO 40
      RHOT = RHOMBI
      GO TO 20
   30 WTFL = WTFL/2.
      NER(2)= NER(2)+1
      WRITE(6,1020) WTFL
      IF(NER(2).EQ.10) STOP
      GO TO 10
   40 VI = RHOVI/RHOMBI
      LAMBDA = RM(MBI)*(VI*SIN(BETAI)+OMEGA*RM(MBI))
C
C  CALCULATE MAXIMUM VALUES FOR RHO*W AT LEADING AND TRAILING EDGE
C
      TWL = 2.*OMEGA*LAMBDA
      AA = (TWL-(OMEGA*RM(MBI))**2)/CPTIP
      TPP = TIP*(1.-AA)
      BB = TGROG*TPP
      TTIP = 1.-BB/CPTIP-AA
      RHOT = RHOIP*TTIP**EXPON
      RHOWMI = RHOT*SQRT(BB)
      AA = (TWL-(OMEGA*RM(MBO2))**2)/CPTIP
      TPP = TIP*(1.-AA)
      BB = TGROG*TPP
      TTIP = 1.-BB/CPTIP-AA
      RHOT = RHOIP*TTIP**EXPON
      RHOWMO = RHOT*SQRT(BB)
C
C  CALCULATE VO AND W-CRITICAL AT BLADE LEADING AND TRAILING EDGE
C
      RHOVO = WTFL/BE(MBO2)/PITCH/COS(BETAO)/RM(MBO2)
      RHOMB2 = RHOIP
      TWLMR = TWL-(OMEGA*RM(MBO2))**2
      LER(1)=1
C     DENSTY CALL NO. 1
      CALL DENSTY(RHOVO,RHOMB2,VO,TWLMR,CPTIP,EXPON,RHOIP,GAM,AR,TIP)
      WCRI = SQRT(TGROG*TIP*(1.-(TWL-(OMEGA*RM(MBI))**2)/CPTIP))
      WCRO = SQRT(TGROG*TIP*(1.-(TWL-(OMEGA*RM(MBO2))**2)/CPTIP))
C
C   CALCULATE BETA CORRECTED TO BOUNDARY A-N AND G-H
C
      TWLMR = TWL-(OMEGA*RM(1))**2
      RHO1 = RHOMBI
      TBI1 = 1.E20
   50 TBIT = (TBI/BE(MBI)*RHO1/RHOMBI+OMEGA*(RM(MBI)**2-RM(1)**2)*RHO1
     1/WTFL*PITCH)*BE(1)
      IF(ABS(TBI1-TBIT).LT..00001) GO TO 60
      TBI1 = TBIT
      RHOVI = WTFL/PITCH*SQRT(1.+TBI1**2)/BE(1)/RM(1)
      LER(1)=2
C     DENSTY CALL NO. 2
      CALL DENSTY (RHOVI,RHO1,AA,TWLMR,CPTIP,EXPON,RHOIP,GAM,AR,TIP)
      GO TO 50
   60 TBI = TBIT
      BTAIN = ATAN(TBI)*57.295779
      TWLMR = TWL-(OMEGA*RM(MM))**2
      RHOMM = RHOMB2
      TBOM = 1.E20
   70 TBOT = (TBO/BE(MBO2)*RHOMM/RHOMB2+OMEGA*(RM(MBO2)**2-RM(MM)**2)*
     1RHOMM/WTFL*PITCH)*BE(MM)
      IF (ABS(TBOM-TBOT).LT..00001) GO TO 80
      TBOM = TBOT
      RHOVO = WTFL/PITCH*SQRT(1.+TBOM**2)/BE(MM)/RM(MM)
      LER(1)=3
C     DENSTY CALL NO. 3
      CALL DENSTY (RHOVO,RHOMM,AA,TWLMR,CPTIP,EXPON,RHOIP,GAM,AR,TIP)
      GO TO 70
   80 TBO = TBOT
      BTAOUT = ATAN(TBO)*57.295779
C
C  CALCULATE TV, ITV, IV, DTDMV, AND BETAV ARRAYS
C
      ITMIN = 0
      ITMAX = NBBI-1
C  TV, ITV, AND DTDMV ON BLADES
      DO 90  IM=MBI,MBO
      LER(2)=1
C     BLCD CALL NO. 1
      CALL BL1(MV(IM),TV(IM,1),DTDMV(IM,1),INF)
      ITV(IM,1)= INT((TV(IM,1)+DTLR)/HT)
      IF (TV(IM,1).GT.-DTLR) ITV(IM,1)=ITV(IM,1)+1
      ITMIN= MIN0(ITMIN,ITV(IM,1))
      LER(2)=2
C     BLCD CALL NO. 2
      CALL BL2(MV(IM),TV(IM,2),DTDMV(IM,2),INF)
      ITV(IM,2)= INT((TV(IM,2)-DTLR)/HT)
      IF (TV(IM,2).LT.DTLR) ITV(IM,2)=ITV(IM,2)-1
   90 ITMAX= MAX0(ITMAX,ITV(IM,2))
      DO 110 IM=MBI2,MBO2
      LER(2)=3
C     BLCD CALL NO. 3
      CALL BL3(MV(IM),TV(IM,3),DTDMV(IM,3),INF)
      ITV(IM,3)= INT((TV(IM,3)+DTLR)/HT)
      IF (TV(IM,3).GT.-DTLR) ITV(IM,3)=ITV(IM,3)+1
      IF (IM.GT.MBO) GO TO 100
      TV(IM,3) = TV(IM,3)+PITCH
  100 ITMIN= MIN0(ITMIN,ITV(IM,3))
      LER(2)=4
C     BLCD CALL NO. 4
      CALL BL4(MV(IM),TV(IM,4),DTDMV(IM,4),INF)
      ITV(IM,4)= INT((TV(IM,4)-DTLR)/HT)
      IF (TV(IM,4).LT.DTLR) ITV(IM,4)=ITV(IM,4)-1
  110 ITMAX= MAX0(ITMAX,ITV(IM,4))
C  ITV AND IV UPSTREAM OF FRONT BLADE
      FIRST = 0
      LAST = NBBI-1
      DO 120 IM=1,MBIM1
      IV(IM) = NBBI*(IM-1)+1
      ITV(IM,1)= FIRST
  120 ITV(IM,2)= LAST
C  ITV BETWEEN FRONT AND REAR BLADES
      IF (MBOP1.GT.MBI2M1) GO TO 140
      LAST= ITV(MBI2,4)
      FIRST= LAST+1-NBBI
      DO 130 IM=MBOP1,MBI2M1
      ITV(IM,3)= FIRST
  130 ITV(IM,4)= LAST
      ITMIN = MIN0(ITMIN,ITV(MBOP1,3))
C  ITV DOWNSTREAM OF REAR BLADE
  140 LAST= ITV(MBO2,4)
      FIRST= LAST+1-NBBI
      DO 150 IM=MBO2P1,MM
      ITV(IM,3)= FIRST
  150 ITV(IM,4)= LAST
      ITMIN = MIN0(ITMIN,ITV(MM,3))
C  FINISH CALCULATING IV ARRAY
      IV(MBI)  = NBBI*MBIM1+1
      MBOT = MIN0(MBO,MBI2M1)
      IF(MBI.GT.MBOT) GO TO 165
      DO 160 IM=MBI,MBOT
  160 IV(IM+1) = IV(IM)+ITV(IM,2)-ITV(IM,1)+1
  165 IF(MBI2.GT.MBO) GO TO 180
      DO 170 IM=MBI2,MBO
  170 IV(IM+1) = IV(IM)+ITV(IM,2)-ITV(IM,3)+ITV(IM,4)-ITV(IM,1)+2-NBBI
  180 DO 190 IM=MBOP1,MM
  190 IV(IM+1) = IV(IM)+ITV(IM,4)-ITV(IM,3)+1
C  BETAV ARRAY
      DO 200 SURF=1,2
      DO 200 IM=MBI,MBO
  200 BETAV(IM,SURF) = ATAN(DTDMV(IM,SURF)*RM(IM))*57.295779
      DO 210 SURF=3,4
      DO 210 IM=MBI2,MBO2
  210 BETAV(IM,SURF) = ATAN(DTDMV(IM,SURF)*RM(IM))*57.295779
      NIP = IV(MM)+NBBI-1
      WRITE(6,1030) VI,RHOWMI,WCRI,BTAIN,VO,RHOWMO,WCRO,BTAOUT
      WRITE(6,1040) PITCH,HT,HM1,HM2,HM3
      WRITE(6,1050) ITMIN,ITMAX,LAMBDA,NIP
      WRITE(6,1060) (SURF,BV(SURF),SURF=1,4)
      IF(BLDAT.LE.0) GO TO 230
      FIRST = MBI
      LAST = MBO
      WRITE (6,1070)
      DO 220 SURF=1,3,2
      I = SURF+1
      WRITE (6,1080) SURF,I,(MV(IM),TV(IM,SURF),DTDMV(IM,SURF),
     1TV(IM,I),DTDMV(IM,I),IM=FIRST,LAST)
      FIRST = MBI2
  220 LAST = MBO2
      WRITE (6,1090) (IM,MV(IM),RM(IM),SAL(IM),BE(IM),DBDM(IM),IM=1,MM)
  230 CONTINUE
C
C  CALCULATE MH AND DTDMH ARRAYS
C
      ITO = ITV(1,1)
      MRTS = 1
      IMS(1) = 1
      MH(1,1) = 0.
      DTDMH(1,1) = 1.E10
      LER(2)=5
C     BLCD AND ROOT (VIA MHORIZ) CALL NO. 5
      CALL MHORIZ(MV,ITV(1,1),BL1,MBI,MBO,ITO,HT,DTLR,0,IMS(1),MH(1,1),
     1DTDMH(1,1),MRTS)
      IF (ITV(MBO,1)-ITV(MBO,2)+NBBI.NE.2) GO TO 240
      IMSL = IMS(1)+1
      MH(IMSL,1) = MV(MBO)
      DTDMH(IMSL,1) = -1.E10
      IMS(1) = IMSL
  240 IMS(2) = 0
      MRTS = 1
      LER(2)=6
C     BLCD AND ROOT (VIA MHORIZ) CALL NO. 6
      CALL MHORIZ(MV,ITV(1,2),BL2,MBI,MBO,ITO,HT,DTLR,1,IMS(2),MH(1,2),
     1DTDMH(1,2),MRTS)
      IMS(3) = 0
      IF (ITV(MBI2,3)-ITV(MBI2,4).NE.2.AND.
     1ITV(MBI2,4)-ITV(MBI2,3).NE.NBBI-2) GO TO 250
      MRTS = 1
      IMS(3) = 1
      MH(1,3) = MV(MBI2)
      DTDMH(1,3) = 1.E10
  250 LER(2)=7
C     BLCD AND ROOT (VIA MHORIZ) CALL NO. 7
      CALL MHORIZ(MV,ITV(1,3),BL3,MBI2,MBO,ITO,HT,DTLR,0,IMS(3),MH(1,3),
     1DTDMH(1,3),MRTS)
      MBOT = MAX0(MBO,MBI2)
      LER(2)=8
C     BLCD AND ROOT (VIA MHORIZ) CALL NO. 8
      CALL MHORIZ(MV,ITV(1,3),BL3,MBOT,MBO2,ITO,HT,DTLR,0,IMS(3),
     1MH(1,3),DTDMH(1,3),MRTS)
      IF (ITV(MBO2,3)-ITV(MBO2,4)+NBBI.NE.2) GO TO 260
      IMSL = IMS(3)+1
      MH(IMSL,3) = MV(MBO2)
      DTDMH(IMSL,3) = -1.E10
      IMS(3) = IMSL
  260 IMS(4) = 0
      IF (ITV(MBI2,3)-ITV(MBI2,4).EQ.2.OR.ITV(MBI2,4)-ITV(MBI2,3).EQ.
     1NBBI-2) MRTS = 1
      LER(2)=9
C     BLCD AND ROOT (VIA MHORIZ) CALL NO. 9
      CALL MHORIZ(MV,ITV(1,4),BL4,MBI2,MBO2,ITO,HT,DTLR,1,IMS(4),
     1MH(1,4),DTDMH(1,4),MRTS)
      I = MAX0(IMS(1),IMS(2),IMS(3),IMS(4))
      IF (I.LE.100) GO TO 270
      WRITE(6,1100) I
      STOP
C
C  CORRECT ITV ARRAY FOR OVERLAPPING PORTION OF BLADE SURFACE 3
C
  270 IF (MBI2.GT.MBO) GO TO 290
      DO 280 IM=MBI2,MBO
  280 ITV(IM,3) = ITV(IM,3)+NBBI
  290 IF(BLDAT.LE.0) GO TO 300
      WRITE (6,1110)  (IM,IV(IM),(ITV(IM,SURF),SURF=1,4),IM=1,MM)
C
C  CALCULATE RMH, BEH, AND BETAH ARRAYS
C
  300 IF(BLDAT.GT.0) WRITE(6,1120)
      DO 320 SURF=1,4
      CALL SPLINT(MR,RMSP,NRSP,MH(1,SURF),IMS(SURF),RMH(1,SURF),AAA)
      CALL SPLINT(MR,BESP,NRSP,MH(1,SURF),IMS(SURF),BEH(1,SURF),AAA)
      IMSS = IMS(SURF)
      IF(IMSS.LT.1) GO TO 320
      DO 310 IHS = 1,IMSS
  310 BETAH(IHS,SURF) = ATAN(DTDMH(IHS,SURF)*RMH(IHS,SURF))*57.295779
      IF (BLDAT.GT.0) WRITE(6,1130) SURF,(MH(IM,SURF),RMH(IM,SURF),
     1BEH(IM,SURF),BETAH(IM,SURF),DTDMH(IM,SURF),IM=1,IMSS)
  320 CONTINUE
      IF (BLDAT.LE.0) GO TO 340
      WRITE (6,1140)
      IT = ITMIN
  330 IF (IT.GT.ITMAX) GO TO 340
      TH = FLOAT(IT)*HT
      WRITE (6,1010) IT,TH
      IT = IT+1
      GO TO 330
  340 IF(NIP.LE.2000) GO TO 350
      WRITE(6,1150)
      STOP
  350 WRITE (6,1000)
      RETURN
 1000 FORMAT (1H1)
 1010 FORMAT (4X,I4,G16.5)
 1020 FORMAT(60HLINPUT WEIGHT FLOW (WTFL) IS TOO LARGE AT BLADE LEADING
     1EDGE/16H WTFL REDUCED TO,G14.6)
 1030 FORMAT (1H1/////24X,10HFREESTREAM,8X,13HMAXIMUM VALUE,
     17X,8HCRITICAL,30X,14HBETA CORRECTED/25X,8HVELOCITY,10X,9HFOR RHO*W
     2,10X,8HVELOCITY,31X,11HTO BOUNDARY/1X,17HLEADING EDGE  B-M,3G18.5,
     312X,12HBOUNDARY A-N,G18.5/1X,17HTRAILING EDGE F-I,3G18.5,12X,
     412HBOUNDARY G-H,G18.5)
 1040 FORMAT (/////5X,28HCALCULATED PROGRAM CONSTANTS//5X,5HPITCH,13X,
     12HHT,13X,3HHM1,13X,3HHM2,13X,3HHM3/1X,5G16.7)
 1050 FORMAT (/5X,5HITMIN,10X,5HITMAX/4X,I5,10X,I5//5X,6HLAMBDA/1X,G16.7
     1///5X,33HNUMBER OF INTERIOR MESH POINTS = ,I5)
 1060 FORMAT (//////5X,23HSURFACE BOUNDARY VALUES//5X,7HSURFACE,7X,2HBV
     1/(5X,I4,4X,F10.5))
 1070 FORMAT (1H1,6X,62HBLADE DATA AT INTERSECTIONS OF VERTICAL MESH LIN
     1ES WITH BLADES)
 1080 FORMAT (1HL,22X,13HBLADE SURFACE,I2,15X,13HBLADE SURFACE,I2/7X,
     1   1HM,14X,2HTV,11X,5HDTDMV,12X,2HTV,11X,5HDTDMV/(5G15.5))
 1090 FORMAT (1H1,13X,44HSTREAM SHEET COORDINATES AND THICKNESS TABLE /
     1   2X,2HIM,7X,1HM,14X,1HR,13X,3HSAL,13X,1HB,12X,5HDB/DM/(1X,I3,
     2   5G15.5))
 1100 FORMAT(34HLONE OF THE MH ARRAYS IS TOO LARGE/7H IT HAS,I5, 8H  POI
     1NTS)
 1110 FORMAT (4H1 IM,9X,8HIV ARRAY,32X,9HITV ARRAY/38X,5HBLADE/37X,7HSUR
     1FACE,3X,1H1,5X,1H2,5X,1H3,5X,1H4/39X,3HNO./(1X,I3,5X,I10,25X,
     24(I4,2X)))
 1120 FORMAT (67H1M COORDINATES OF INTERSECTIONS OF HORIZONTAL MESH LINE
     1S WITH BLADE)
 1130 FORMAT (25HLMH ARRAY - BLADE SURFACE,I2//15X,2HMH,19X,3HRMH,19X,
     1   3HBEH,18X,5HBETAH,17X,5HDTDMH/(5G22.4))
 1140 FORMAT (43H1THETA COORDINATES OF HORIZONTAL MESH LINES//6X,2HIT,
     15X,5HTHETA)
 1150 FORMAT(48HLTHE NUMBER OF INTERIOR MESH POINTS EXCEEDS 2000)
      END
$IBFTC MHORIZ  DEBUG
      SUBROUTINE MHORIZ(MV,ITV,BL,MBI,MBO,ITO,HT,DTLR,KODE,J,MH,DTDMH,
     1MRTS)
C
C  MHORIZ CALCULATES M COORDINATES OF INTERSECTIONS OF ALL HORIZONTAL
C  MESH LINES WITH A BLADE SURFACE
C  KODE = 0 FOR UPPER BLADE SURFACE
C  KODE = 1 FOR LOWER BLADE SURFACE
C
      COMMON SRW,ITER,IEND,LER(2),NER(2)
      DIMENSION MV(100),ITV(100),MH(100),DTDMH(100)
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      REAL MVIM
      EXTERNAL BL
      IF (MBI.GE.MBO) RETURN
      IM= MBI
   10 ITIND= 0
   20 IF (ITV(IM+1)-ITV(IM)-ITIND) 30,40,50
   30 J= J+1
      TI= FLOAT(ITV(IM+1)-ITO-ITIND+KODE)*HT
      ITIND= ITIND-1
      MVIM = MV(IM)
      IF (MRTS.EQ.1) MVIM = MVIM+(MV(IM+1)-MVIM)/1000.
      CALL ROOT (MVIM,MV(IM+1),TI,BL,DTLR,MH(J),DTDMH(J))
      GO TO 20
   40 IM= IM+1
      MRTS = 0
      IF (IM.EQ.MBO) RETURN
      GO TO 10
   50 J= J+1
      TI= FLOAT(ITV(IM)-ITO+ITIND+KODE)*HT
      ITIND= ITIND+1
      MVIM = MV(IM)
      IF (MRTS.EQ.1) MVIM = MVIM+(MV(IM+1)-MVIM)/1000.
      CALL ROOT (MVIM,MV(IM+1),TI,BL,DTLR,MH(J),DTDMH(J))
      GO TO 20
      END
$IBFTC COEF    DEBUG
      SUBROUTINE COEF
C
C  COEF CALCULATES FINITE DIFFERENCE COEFFICIENTS, A, AND CONSTANTS, K,
C  AT ALL UNKNOWN MESH POINTS FOR THE ENTIRE REGION
C
      COMMON SRW,ITER,IEND,LER(2),NER(2)
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /HRBAAK/ H(4),R(4),B(4),KAK(4),KA(4),IH(4),RZ,BZ
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
C  INITIALIZE ARRAYS
      ITER = ITER+1
      IH(1) = 1
      DO 10 I=2,4
   10 IH(I) = 0
      IF(ITV(MBI2,3)-ITV(MBI2,4).EQ.2) IH(3) = 1
      IF(ITV(MBI2,4)-ITV(MBI2,3).EQ.NBBI-2.AND.MBI2.NE.MBOP1) IH(3) = 1
C  INCOMPRESSIBLE CASE
      IF(GAM.NE.1.5.OR.AR.NE.1000..OR.TIP.NE.1.E6) GO TO 20
      IEND = 1
      GO TO 40
C  ADJUSTMENT OF PRINTING CONTROL VARIABLES
   20 IF(ITER.NE.1.AND.ITER.NE.2) GO TO 30
      AANDK = AANDK-1
      ERSOR = ERSOR-1
      STRFN = STRFN-1
      SLCRD = SLCRD-1
      INTVL = INTVL-1
      SURVL = SURVL-1
   30 IF(IEND.NE.0) GO TO 40
      AANDK = AANDK+2
      ERSOR = ERSOR+2
      STRFN = STRFN+2
      SLCRD = SLCRD+2
      INTVL = INTVL+2
      SURVL = SURVL+2
C
C  FIRST VERTICAL MESH LINE
C
   40 DO 50 IP=1,NBBI
      A(IP,1) = 0.
      A(IP,2) = 0.
      A(IP,3) = 0.
      A(IP,4) = 1.
   50 K(IP) = HM1*TBI/PITCH/RM(1)
C
C  UPSTREAM OF BLADES, EXCEPT FOR FIRST VERTICAL MESH LINE
C
      IF(2.GT.MBIM1) GO TO 70
      DO 60 IM=2,MBIM1
   60 CALL COEFP(IM,1,2)
C
C  BETWEEN FRONT BLADES
C
   70 MBOT = MIN0(MBO,MBI2M1)
      IF(MBI.GT.MBOT) GO TO 90
      DO 80 IM=MBI,MBOT
   80 CALL COEFBB(IM,1,2,1)
   90 IF(MBI2.GT.MBO) GO TO 110
C
C  OVERLAPPING CASE
C
      DO 100 IM=MBI2,MBO
      CALL COEFBB(IM,1,4,1)
  100 CALL COEFBB(IM,3,2,4)
      GO TO 130
C
C   NON - OVERLAPPING CASE
C
  110 IF(MBOP1.GT.MBI2M1) GO TO 130
      DO 120 IM=MBOP1,MBI2M1
  120 CALL COEFP(IM,3,4)
C
C  BETWEEN REAR BLADES
C
  130 MBIT = MAX0(MBI2,MBOP1)
      IF(MBIT.GT.MBO2) GO TO 150
      DO 140 IM=MBIT,MBO2
  140 CALL COEFBB(IM,3,4,3)
C
C  DOWNSTREAM OF BLADES EXCEPT FOR FINAL MESH LINE
C
  150 IF(MBO2P1.GT.MMM1) GO TO 170
      DO 160 IM=MBO2P1,MMM1
  160 CALL COEFP(IM,3,4)
C
C  FINAL VERTICAL MESH LINE
C
  170 IVMM = IV(MM)
      DO 180 IP=IVMM,NIP
      A(IP,1) = 0.
      A(IP,2) = 0.
      A(IP,3) = 1.
      A(IP,4) = 0.
  180 K(IP) = -HM3*TBO/PITCH/RM(MM)
C
C  TAKE CARE OF POINTS ADJACENT TO B, AND CASES WHEN POINTS J,C,E, OR F
C  ARE GRID POINTS
C
C  POINT B
      IP = IV(MBIM1)
      A(IP,4) = 0.
C  POINT J
      IF(ITV(MBI2,3)-ITV(MBI2,4).NE.2) GO TO 190
      IT = ITV(MBI2,4)+1
      IP = IPF(MBI2M1,IT)
      K(IP) = K(IP)+A(IP,4)*BV(4)
      A(IP,4) = 0.
C  POINT C
  190 IF(ITV(MBO,1)-ITV(MBO,2)+NBBI.NE.2) GO TO 200
      IT = ITV(MBO,1)-1
      IP = IPF(MBOP1,IT)
      A(IP,3) = 0.
C  POINT E
  200 IF(ITV(MBI2,4)-ITV(MBI2,3).NE.NBBI-2.OR.MBI2.EQ.MBOP1) GO TO 210
      IP = IV(MBI2M1)
      K(IP) = K(IP)+A(IP,4)*BV(3)
      A(IP,4) = 0.
C  POINT F
  210 IF ((ITV(MBO2,3)-ITV(MBO2,4)+NBBI.NE.2).AND.(ITV(MBO2,3)
     1-ITV(MBO2,4).NE.2)) GO TO 220
      IP = IV(MBO2P1)
      K(IP) = K(IP)+A(IP,3)*BV(3)
      A(IP,3) = 0.
C
C   LINE K-L AND LINE TO RIGHT OF C-D
C
  220 FIRST = MAX0(ITV(MBOP1,3)+NBBI,ITV(MBO,3))
      IPKL = IPF(MBO,FIRST)
      IPL = IPKL+ITV(MBO,2)-FIRST
      IPCD = IPF(MBOP1,FIRST-NBBI)
  230 IF(IPKL.GT.IPL) RETURN
      K(IPKL) = K(IPKL)+A(IPKL,4)
      K(IPCD) = K(IPCD)-A(IPCD,3)
      IPKL = IPKL+1
      IPCD = IPCD+1
      GO TO 230
      END
$IBFTC COEFBB  DEBUG
      SUBROUTINE COEFBB(IM,UPPER,LOWER,UPPRBV)
C
C  COEFBB CALCULATES FINITE DIFFERENCE COEFFICIENTS, A, AND CONSTANTS, K
C  ALONG ALL VERTICAL MESH LINES WHICH INTERSECT BLADES
C
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /HRBAAK/ H(4),R(4),B(4),KAK(4),KA(4),IH(4),RZ,BZ
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      IF(ITV(IM,UPPER).GT.ITV(IM,LOWER)) RETURN
      ITVU = ITV(IM,UPPER)
      ITVL = ITV(IM,LOWER)
      IT   = ITVU - 1
      IPU = IPF(IM,ITVU)
      IPL = IPU+ITVL-ITVU
      DO 90 IP=IPU,IPL
      IT = IT+1
      CALL HRB(IM,IT,IP)
      DO 10 I=1,4
      KAK(I) = 0.
   10 KA(I) = 0
C  FIX HRB VALUES FOR LINES C-D AND K-L
      IF(IM.NE.MBOP1) GO TO 20
      IF(IT.GE.ITV(IM-1,1)) GO TO 20
      IP3 = IPF(IM-1,IT+NBBI)
      R(3) = RHO(IP3)
   20 IF(IM+1.NE.MBOP1) GO TO 60
      IF(MBI2P1-MBOP1) 30,30,40
   30 IF(IT-ITV(IM,3)) 60,50,50
   40 IF(IT.LE.ITV(IM+1,4)) GO TO 60
   50 IP4 = IPF(IM+1,IT-NBBI)
      R(4) = RHO(IP4)
C  FIX HRB VALUES FOR CASES WHERE MESH LINES INTERSECT BLADES
   60 IF (IT .EQ. ITV(IM,UPPER)) CALL BDRY12(1,IM,IT,UPPER,UPPRBV)
      IF (IT .EQ. ITV(IM,LOWER)) CALL BDRY12(2,IM,IT,LOWER,LOWER)
      ITVM1 = ITV(IM-1,UPPER)
      ITVP1 = ITV(IM+1,UPPER)
      IF (IM.EQ.MBO.AND.UPPER.EQ.3) ITVP1 = ITVP1+NBBI
      IF (IM.EQ.MBOP1) ITVM1 = ITVM1-NBBI
      IF (IT.LT.ITVM1) CALL BDRY34(3,IM,UPPER,UPPRBV)
      IF (IT.LT.ITVP1) CALL BDRY34(4,IM,UPPER,UPPRBV)
      IF(IM.EQ.MBI2.AND.LOWER.EQ.4) GO TO 70
      IF(IM.EQ.MBOP1.AND.LOWER.EQ.4.AND.MBI2.GT.MBO) GO TO 70
      IF (IT.GT.ITV(IM-1,LOWER)) CALL BDRY34(3,IM,LOWER,LOWER)
   70 IF(IM.EQ.MBO.AND.LOWER.EQ.2) GO TO 80
      IF (IT.GT.ITV(IM+1,LOWER)) CALL BDRY34(4,IM,LOWER,LOWER)
C  COMPUTE A AND K COEFFICIENTS
   80 CALL AAK(IM,IP)
      DO 90 I=1,4
      K(IP) = K(IP)+KAK(I)*A(IP,I)
   90 IF(KA(I).EQ.1) A(IP,I) = 0.
      RETURN
C
C  COEFP CALCULATES FINITE DIFFERENCE COEFFICIENTS, A, AND CONSTANTS, K,
C  ALONG ALL VERTICAL MESH LINES WHICH DO NOT INTERSECT BLADES
C
      ENTRY COEFP(IM,UPPER,LOWER)
      ITVU = ITV(IM,UPPER)
      ITVL = ITV(IM,LOWER)
      IT = ITVU-1
      IPU = IV(IM)
      IPL = IV(IM+1)-1
      DO 100 IP=IPU,IPL
      IT = IT+1
      CALL HRB(IM,IT,IP)
      IF (IT.EQ.ITVU) R(1) = RHO(IPL)
      IF (IT.EQ.ITVL) R(2) = RHO(IPU)
      IF(IM.NE.MBOP1) GO TO 100
      IF(IT.GE.ITV(IM-1,1)) GO TO 100
      IP3 = IPF(IM-1,IT+NBBI)
      R(3) = RHO(IP3)
  100 CALL AAK(IM,IP)
      K(IPL) = K(IPL)+A(IPL,2)
      K(IPU) = K(IPU)-A(IPU,1)
      RETURN
      END
$IBFTC HRB     DEBUG
      SUBROUTINE HRB(IM,IT,IP)
C
C  HRB CALCULATES MESH SPACING, H, DENSITIES, RZ AND R, AT GIVEN AND
C  ADJACENT POINTS, AND STREAM SHEET THICKNESSES, BZ AND B, AT GIVEN
C  AND ADJACENT POINTS
C
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /HRBAAK/ H(4),R(4),B(4),KAK(4),KA(4),IH(4),RZ,BZ
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      H(1)= HT*RM(IM)
      H(2)= HT*RM(IM)
      H(3)= MV(IM)-MV(IM-1)
      H(4)= MV(IM+1)-MV(IM)
      RZ = RHO(IP)
      IP3 = IPF(IM-1,IT)
      IP4 = IPF(IM+1,IT)
      R(1) = RHO(IP-1)
      R(2) = RHO(IP+1)
      R(3) = RHO(IP3)
      R(4) = RHO(IP4)
      BZ= BE(IM)
      B(3)= BE(IM-1)
      B(4)= BE(IM+1)
      RETURN
      END
$IBFTC AAK     DEBUG
      SUBROUTINE AAK(IM,IP)
C
C  AAK CALCULATES FINITE DIFFERENCE COEFFICIENTS, A, AND CONSTANT, K,
C  AT A SINGLE MESH POINT
C
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /HRBAAK/ H(4),R(4),B(4),KAK(4),KA(4),IH(4),RZ,BZ
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      A12= 2./H(1)/H(2)
      A34= 2./H(3)/H(4)
      AZ= A12+A34
      B12= (R(2)-R(1))/RZ/(H(1)+H(2))
      B34= (B(4)*R(4)-B(3)*R(3))/BZ/RZ/(H(3)+H(4))-SAL(IM)/RM(IM)
      A(IP,1) = (2./H(1)+B12)/AZ/(H(1)+H(2))
      A(IP,2) = A12/AZ-A(IP,1)
      A(IP,3) = (2./H(3)+B34)/AZ/(H(3)+H(4))
      A(IP,4) = A34/AZ-A(IP,3)
      K(IP) = -TWW*BZ*RZ*SAL(IM)/AZ
      RETURN
      END
$IBFTC BDRY12  DEBUG
      SUBROUTINE BDRY12(I,IM,IT,SURF,SURFBV)
C
C  BDRY12 CORRECTS VALUES COMPUTED BY HRB WHEN A VERTICAL MESH LINE
C  INTERSECTS A BLADE
C
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /RHOS/ RHOHB(100,4),RHOVB(100,4)
      COMMON /HRBAAK/ H(4),R(4),B(4),KAK(4),KA(4),IH(4),RZ,BZ
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      H(I) = ABS(FLOAT(IT)*HT-TV(IM,SURF))*RM(IM)
      R(I)= RHOVB(IM,SURF)
      KAK(I)=BV(SURFBV)
      KA(I)=1
      RETURN
      END
$IBFTC BDRY34  DEBUG
      SUBROUTINE BDRY34(I,IM,SURF,SURFBV)
C
C  BDRY34 CORRECTS VALUES COMPUTED BY HRB WHEN A HORIZONTAL MESH LINE
C  INTERSECTS A BLADE
C
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /RHOS/ RHOHB(100,4),RHOVB(100,4)
      COMMON /HRBAAK/ H(4),R(4),B(4),KAK(4),KA(4),IH(4),RZ,BZ
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      IH(SURF)=IH(SURF)+1
      IHS=IH(SURF)
      H(I)=ABS(MV(IM)-MH(IHS,SURF))
      R(I)=RHOHB(IHS,SURF)
      B(I)=BEH(IHS,SURF)
      KAK(I)=BV(SURFBV)
      KA(I)=1
      RETURN
      END
$IBFTC SOR     DEBUG
      SUBROUTINE SOR
C
C  SOR SOLVES THE SET OF SIMULTANEOUS EQUATIONS FOR THE STREAM FUNCTION
C  USING THE METHOD OF SUCCESSIVE OVER-RELAXATION
C
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      AATEMP = AANDK
      IF (ORF.GE.2.) ORF=0.
      IF (ORF.GT.1.) GO TO 20
      ORF = 1.
      ORFOPT= 2.
   10 ORFTEM=ORFOPT
      LMAX = 0.
   20 IF (AATEMP.GT.0) WRITE(6,1010)
      ERROR = 0.
C
C  SOLVE MATRIX EQUATION BY SOR, OR CALCULATE OPTIMUM OVERRELAXATION
C  FACTOR
C
      IP = 0
      DO 100 IM=1,MM
      IPU = IV(IM)
      IPL = IV(IM+1)-1
      IT = ITV(IM,1)
      IF (IM.GT.MBO) IT=ITV(IM,3)
      IF (AATEMP.GT.0) WRITE(6,1020) IM,IT
      DO 100 IP=IPU,IPL
      IF(IT.GT.ITV(IM,4).AND.IT.LE.ITV(IM,3)) IT = IT+ITV(IM,3)
     1-ITV(IM,4)-1
      IP1 = IP-1
      IP2 = IP+1
C  CORRECT IP1 AND IP2 ALONG PERIODIC BOUNDARIES
      IF(IM.GE.MBI.AND.IM.LE.MBO2.AND.(IM.LE.MBO.OR.IM.GE.MBI2))GO TO 30
      IF(IT.EQ.ITV(IM,1).OR.IT.EQ.ITV(IM,3)) IP1 = IP1+NBBI
      IF(IT.EQ.ITV(IM,2).OR.IT.EQ.ITV(IM,4)) IP2 = IP2-NBBI
   30 IT3 = IT
      IT4 = IT
C  CORRECT IT3 AND IT4 ALONG LINES C-D AND K-L
      IF(IM.NE.MBOP1) GO TO 40
      IF(IT.LT.ITV(IM-1,1)) IT3 = IT+NBBI
   40 IF(IM.NE.MBO) GO TO 70
      IF(MBI2-MBO) 50,50,60
   50 IF(IT.GE.ITV(IM,3)) IT4 = IT-NBBI
      GO TO 70
   60 IF (IT.GT.ITV(IM+1,4)) IT4 = IT-NBBI
   70 IP3 = IPF(IM-1,IT3)
      IP4 = IPF(IM+1,IT4)
      IF (ORF.GT.1.) GO TO 80
C  CALCULATE NEW ESTIMATE FOR LMAX
      UNEW = A(IP,1)*U(IP1)+A(IP,2)*U(IP2)+A(IP,3)*U(IP3)+A(IP,4)*U(IP4)
      IF (UNEW.LT.1.E-25) U(IP) = 0.
      IF (U(IP).EQ.0.) GO TO 90
      RATIO = UNEW/U(IP)
      LMAX= AMAX1(RATIO,LMAX)
      U(IP) = UNEW
      GO TO 90
C  CALCULATE NEW ESTIMATE FOR STREAM FUNCTION BY SOR
   80 CHANGE = ORF*(K(IP)-U(IP)+A(IP,1)*U(IP1)+A(IP,2)*U(IP2)+A(IP,3)*
     1U(IP3)+A(IP,4)*U(IP4))
      ERROR= AMAX1(ERROR,ABS(CHANGE))
      U(IP) = U(IP)+CHANGE
   90 IF(AATEMP.LE.0) GO TO 100
      WRITE (6,1030) IT,IP,IP1,IP2,IP3,IP4,(A(IP,I),I=1,4),K(IP)
  100 IT = IT+1
      AATEMP = 0
      IF(ORF.GT.1.) GO TO 110
      ORFOPT= 2./(1.+SQRT(ABS(1.-LMAX)))
      WRITE(6,1040) ORFOPT
      IF(ORFTEM-ORFOPT.GT..00001.OR.ORFOPT.GT.1.999) GO TO 10
      WRITE (6,1000)
      ORF = ORFOPT
      GO TO 20
  110 IF(ERSOR.GT.0) WRITE(6,1050) ERROR
      IF(ERROR.GT..000001) GO TO 20
      IF(STRFN.LE.0) RETURN
C
C  PRINT STREAM FUNCTION VALUES FOR THIS ITERATION
C
      WRITE (6,1060)
      MBIT = MIN0(MBO,MBI2M1)
      DO 120 IM =1,MBIT
      IPU = IV(IM)
      IPL = IV(IM+1)-1
      ITVU = ITV(IM,1)
      WRITE (6,1020) IM,ITVU
  120 WRITE (6,1070) (U(IP),IP=IPU,IPL)
      IF(MBI2.GT.MBO) GO TO 140
      DO 130 IM=MBI2,MBO
      IPU = IV(IM)
      IPL = IV(IM)+ITV(IM,4)-ITV(IM,1)
      ITVU = ITV(IM,1)
      WRITE (6,1020) IM,ITVU
      WRITE (6,1070) (U(IP),IP=IPU,IPL)
      IPU = IPL+1
      IPL = IV(IM+1)-1
      IF(IPU.GT.IPL) GO TO 130
      ITVU = ITV(IM,3)
      WRITE (6,1020) IM,ITVU
      WRITE (6,1070) (U(IP),IP=IPU,IPL)
  130 CONTINUE
  140 DO 150 IM=MBOP1,MM
      IPU = IV(IM)
      IPL = IV(IM+1)-1
      ITVU = ITV(IM,3)
      WRITE (6,1020) IM,ITVU
  150 WRITE (6,1070) (U(IP),IP=IPU,IPL)
      RETURN
 1000 FORMAT (1H1)
 1010 FORMAT (82H1  IT   IP    IP1   IP2   IP3   IP4    A(1)      A(2)
     1    A(3)      A(4)        K)
 1020 FORMAT(5H IM =,I4,6X,6HIT1 = ,I4)
 1030 FORMAT(1X,I4,5I6,5F10.5)
 1040 FORMAT(24H ESTIMATED OPTIMUM ORF =,F9.6)
 1050 FORMAT(8H ERROR =,F11.8)
 1060 FORMAT(1H1,10X,22HSTREAM FUNCTION VALUES)
 1070 FORMAT (2X,10F13.8)
      END
$IBFTC SLAX    DEBUG
      SUBROUTINE SLAX
C
C  SLAX CALLS SUBROUTINES TO CALCULATE RHO*W-SUB-M THROUGHOUT THE REGION
C  AND ON THE BLADE SURFACES, AND TO CALCULATE AND PLOT THE
C  STREAMLINE LOCATIONS
C
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /SLA/ TSL(800),UINT(8)
      DIMENSION MSL(100),KKK(18),P(4)
      DIMENSION W(2000),RWM(2000),BETA(2000),WMB(100,4),WTB(100,4),
     1XDOWN(800),YACROS(800)
      EQUIVALENCE (A(1,1),W(1)),(A(1,2),RWM(1)),(A(1,3),BETA(1)),
     1(A(1,4),WMB(1)),(A(401,4),WTB(1)),(A(801,4),XDOWN(1)),
     2(K(1),YACROS(1))
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      DATA (KKK(J),J=4,18,2)/8*1H*/
C
C  CALL SLAVP AND SLAVBB THROUGHOUT THE REGION
C
      ITVU= ITV(1,1)
      ITVL= ITV(1,2)
      DO 10 IM=1,MBIM1
   10 CALL SLAVP(IM,ITVU,ITVL)
      MBOT= MIN0(MBO,MBI2M1)
      IF (MBI.GT.MBOT) GO TO 30
      DO 20 IM=MBI,MBOT
      I= 0
   20 CALL SLAVBB(IM,1,2,1,I)
   30 IF (MBOP1.GT.MBI2M1) GO TO 50
      ITVU= ITV(MBOP1,3)
      ITVL= ITV(MBOP1,4)
      DO 40 IM=MBOP1,MBI2M1
   40 CALL SLAVP(IM,ITVU,ITVL)
   50 IF (MBI2.GT.MBO) GO TO 70
      DO 60 IM=MBI2,MBO
      I= 0
      CALL SLAVBB(IM,1,4,1,I)
   60 CALL SLAVBB(IM,3,2,4,I)
   70 MBOT= MAX0(MBOP1,MBI2)
      IF (MBOT.GT.MBO2) GO TO 90
      DO 80 IM= MBOT,MBO2
      I= 0
   80 CALL SLAVBB(IM,3,4,3,I)
   90 ITVU= ITV(MBO2P1,3)
      ITVL= ITV(MBO2P1,4)
      DO 100 IM=MBO2P1,MM
  100 CALL SLAVP(IM,ITVU,ITVL)
C
C  PLOT STREAMLINES
C
      IF (SLCRD.LE.0) RETURN
      DO 110 IM=1,MM
  110 MSL(IM) = MV(IM)
      KKK(1) = 7
      KKK(2) = 8
      KKK(3) = MM
      P(1) = 1.
      P(3) = 0.
      P(4) = 0.
      WRITE(6,1000)
      CALL PLOTMY(MSL,TSL,KKK,P)
      WRITE(6,1010)
      RETURN
 1000 FORMAT (2HPT,50X,16HSTREAMLINE PLOTS  )
 1010 FORMAT (2HPL,40X,70HSTREAMLINES ARE PLOTTED WITH THETA ACROSS THE
     1PAGE AND M DOWN THE PAGE)
      END
$IBFTC SLAV    DEBUG
      SUBROUTINE SLAV
C
C  SLAV CALCULATES RHO*W-SUB-M THROUGHOUT THE REGION AND ON THE BLADE
C  SURFACES, AND CALCULATES THE STREAMLINE LOCATIONS
C
      COMMON SRW,ITER,IEND,LER(2),NER(2)
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /SLA/ TSL(800),UINT(8)
      DIMENSION TSP(50),USP(50),DUDT(50),TINT(8)
      DIMENSION W(2000),RWM(2000),BETA(2000),WMB(100,4),WTB(100,4),
     1XDOWN(800),YACROS(800)
      EQUIVALENCE (A(1,1),W(1)),(A(1,2),RWM(1)),(A(1,3),BETA(1)),
     1(A(1,4),WMB(1)),(A(401,4),WTB(1)),(A(801,4),XDOWN(1)),
     2(K(1),YACROS(1))
C
C  SLAVP CALCULATES ALONG VERTICAL MESH LINES WHICH DO NOT
C  INTERSECT BLADES
C
      ENTRY SLAVP(IM,ITVU,ITVL)
      LOC= 0
      I1 = 0
      NI= 8
      NSP= ITVL-ITVU+2
      IP = IV(IM)-1
      DO 10 IT=1,NSP
      IP = IP+1
      TSP(IT) = FLOAT(IT+ITVU-1)*HT
   10 USP(IT)= U(IP)
      USP(NSP) = USP(1)+1.
      IP = IV(IM)
      INTU = INT(U(IP)*5.)
      IF (U(IP).GT.0.) INTU=INTU+1
      DO 20 J=1,5
      UINT(J) = FLOAT(INTU)/5.
   20 INTU = INTU+1
      UINT(6)= BV(4)
   30 IF (UINT(6).GE.U(IP)) GO TO 40
      UINT(6)= UINT(6)+1.
      GO TO 30
   40 IF (UINT(6).LE.U(IP)+1.) GO TO 50
      UINT(6)= UINT(6)-1.
      GO TO 40
   50 UINT(7)= UINT(1)
      UINT(8)= UINT(1)
      GO TO 100
C
C  SLAVBB CALCULATES ALONG VERTICAL MESHLINES WHICH INTERSECT BLADES
C
      ENTRY SLAVBB(IM,UPPER,LOWER,UPPRBV,I)
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      LOC= 1
      ITVUP1 = ITV(IM,UPPER)
      ITVLM1 = ITV(IM,LOWER)
      ITVU = ITVUP1-1
      ITVL = ITVLM1+1
      NSP = ITVL-ITVU+1
      TSP(1) = TV(IM,UPPER)
      TSP(NSP) = TV(IM,LOWER)
      USP(1) = BV(UPPRBV)
      USP(NSP) = BV(LOWER)
      IP = IPF(IM,ITVUP1)-1
      NSPM1 = NSP-1
      IF(2.GT.NSPM1) GO TO 70
      DO 60 IT=2,NSPM1
      IP = IP+1
      TSP(IT) = FLOAT(IT+ITVU-1)*HT
   60 USP(IT) = U(IP)
   70 I1= I
      I= I1+1
      UINT(I)= BV(UPPRBV)
      INTU = INT(UINT(I)*5.)
      IF (UINT(I).GE.0.) INTU=INTU+1
   80 I = I+1
      UINT(I) = FLOAT(INTU)/5.
      INTU = INTU+1
      IF (UINT(I).LT.BV(LOWER)) GO TO 80
      UINT(I)= BV(LOWER)
      IF(LOWER-UPPER.NE.1) GO TO 90
      I = 7
      UINT(I) = BV(4)
   90 UINT(8) = BV(LOWER)
      NI= I-I1
      IF (NI.EQ.7) NI = 8
C
C  FOR BOTH SLAVP AND SLAVBB, CALCULATE RHO*W-SUB-M IN THE REGION, AND
C  RHO*W AT VERTICAL MESH LINE INTERSECTIONS ON THE BLADE SURFACES
C
  100 CALL SPLINE(TSP,USP,NSP,DUDT,AAA)
      FIRST= (1-LOC)*ITVU+LOC*ITVUP1
      LAST = (1-LOC)*ITVL+LOC*ITVLM1
      IF(FIRST.GT.LAST) GO TO 120
      IT = LOC
      IPU = IPF(IM,FIRST)
      IPL = IPF(IM,LAST)
      DO 110 IP=IPU,IPL
      IT = IT+1
  110 RWM(IP)  = DUDT(IT)*WTFL/BE(IM)/RM(IM)
  120 IF (LOC.EQ.0) GO TO 130
      WMB(IM,UPPER) = DUDT(  1)*WTFL/BE(IM)/RM(IM)
      WMB(IM,LOWER) = DUDT(NSP)*WTFL/BE(IM)/RM(IM)
      RMDTU2 = (RM(IM)*DTDMV(IM,UPPER))**2
      RMDTL2 = (RM(IM)*DTDMV(IM,LOWER))**2
      IF (RMDTU2.GT.10000.) WMB(IM,UPPER) = 0.
      IF (RMDTL2.GT.10000.) WMB(IM,LOWER) = 0.
      WMB(IM,UPPER) = ABS(WMB(IM,UPPER))*SQRT(1.+RMDTU2)
      WMB(IM,LOWER) = ABS(WMB(IM,LOWER))*SQRT(1.+RMDTL2)
  130 IF (SLCRD.LE.0) RETURN
      I1 = I1+1
      CALL SPLINT(USP,TSP,NSP,UINT(I1),NI,TINT(I1),AAA)
      I2 = NI+I1-1
      DO 140 J=I1,I2
      L= (J-1)*MM+IM
  140 TSL(L)= TINT(J)
      IF (UPPER.EQ.1.AND.LOWER.EQ.4) RETURN
      IF (IM.EQ.1) WRITE(6,1000)
      WRITE(6,1010) MV(IM),(UINT(J),TINT(J),J=1,8)
      RETURN
 1000 FORMAT(1H1////30X,22HSTREAMLINE COORDINATES////5X,12HM COORDINATE,
     13(7X,10HSTREAM FN.,10X,5HTHETA,4X)//)
 1010 FORMAT(1X,7G18.7/(19X,6G18.7))
      END
$IBFTC TANG    DEBUG
      SUBROUTINE TANG
C
C  TANG CALCULATES RHO*W-SUB-THETA AND THEN RHO*W THROUGHOUT THE REGION
C  AND ON THE BLADE SURFACES, AND CALCULATES THE VELOCITY ANGLE, BETA,
C  THROUGHOUT THE REGION
C
      COMMON SRW,ITER,IEND,LER(2),NER(2)
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      DIMENSION SPM(100),USP(100),DUDM(100)
      DIMENSION W(2000),RWM(2000),BETA(2000),WMB(100,4),WTB(100,4),
     1XDOWN(800),YACROS(800)
      EQUIVALENCE (A(1,1),W(1)),(A(1,2),RWM(1)),(A(1,3),BETA(1)),
     1(A(1,4),WMB(1)),(A(401,4),WTB(1)),(A(801,4),XDOWN(1)),
     2(K(1),YACROS(1))
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      LOGICAL ADDL,ADD
      EXTERNAL BL1,BL2,BL3,BL4
C
C  PERFORM CALCULATIONS ALONG ONE HORIZONTAL LINE AT A TIME
C
      IT = ITMIN
   10 IF (IT.GT.ITMAX) RETURN
      ADDL = .FALSE.
      ADD  = .FALSE.
      ITI = IT
      S1 = 0
C
C  ON THE GIVEN HORIZONTAL MESH LINE, FIND A FIRST POINT IN THE REGION
C
      IF(IT.GE.0.AND.IT.LT.NBBI) GO TO 60
      IM = MBIM1
   20 IM = IM+1
      ITV1 = ITV(MBO,3)-NBBI
      IF(MBI2.GT.MBO) ITV1 = ITV(MBOP1,3)
      IF(IM.EQ.MBO.AND.IT.GE.ITV1.AND.IT.LE.ITV(MBO,2)-NBBI) GO TO 200
      IF(IM.GT.MBO2P1) GO TO 200
      DO 40 SURF=1,3,2
      IF (IM.GT.MBO.AND.SURF.EQ.1) GO TO 40
      ITVIM1 = ITV(IM-1,SURF)
      IF(SURF.EQ.3.AND.IM.EQ.MBOP1) ITVIM1 = ITVIM1-NBBI
      IF(ITI.GE.ITV(IM,SURF).AND.ITI.LT.ITVIM1) GO TO 70
   40 CONTINUE
      SURF = 1
      IF(IM.EQ.MBOP1.AND.IT.EQ.ITV(MBO,1)-1.AND.ITV(MBO,1)-ITV(MBO,2)
     1+NBBI.EQ.2) GO TO 70
      DO 50 SURF=2,4,2
      IF (IM.LE.MBI2.AND.SURF.EQ.4) GO TO 50
      IF (IT.LE.ITV(IM,SURF).AND.IT.GT.ITV(IM-1,SURF)) GO TO 70
   50 CONTINUE
      GO TO 20
C
C  FIRST POINT IS ON BOUNDARY A-N
C
   60 IM1= 1
      IM = 1
      SPM(1) = MV(1)
      USP (1) = U(IT+1)
      GO TO 80
C
C  FIRST POINT IS ON A BLADE SURFACE
C
   70 S1 = SURF
      ITI = IT
      ADD = .FALSE.
      IF (ADDL.AND.S1.EQ.3) ADD = .TRUE.
      IF (ADD) ITI=IT-NBBI
      IM1 = IM-1
      IM2 = IM
      TH = FLOAT(ITI)*HT
      IF (S1.EQ.3.AND.IM1.LT.MBO) TH = TH-FLOAT(NBBI)*HT
      MVIM1 = MV(IM1)
      IF ((IM.EQ.MBIP1.AND.(SURF.EQ.1.OR.SURF.EQ.2)).OR.(IM.EQ.MBI2P1
     1.AND.(SURF.EQ.3.OR.SURF.EQ.4))) MVIM1=MVIM1+(MV(IM2)-MVIM1)/1000.
      LER(2)=10
C     BLCD (VIA ROOT) CALL NO. 10
      IF (S1.EQ.1.AND.IM1.NE.MBO) CALL ROOT(MVIM1,MV(IM2),TH,BL1,DTLR,
     1ANS,AAA)
      LER(2)=11
C     BLCD (VIA ROOT) CALL NO. 11
      IF (S1.EQ.3.AND.IM1.NE.MBO2) CALL ROOT(MVIM1,MV(IM2),TH,BL3,DTLR,
     1ANS,AAA)
      LER(2)=12
C     BLCD (VIA ROOT) CALL NO. 12
      IF (S1.EQ.2) CALL ROOT(MVIM1,MV(IM2),TH,BL2,DTLR,ANS,AAA)
      LER(2)=13
C     BLCD (VIA ROOT) CALL NO. 13
      IF(S1.EQ.4) CALL ROOT(MVIM1,MV(IM2),TH,BL4,DTLR,ANS,AAA)
      IF(S1.EQ.1.AND.IM1.EQ.MBO) ANS = MV(MBO)
      IF(S1.EQ.3.AND.IM1.EQ.MBO2) ANS = MV(MBO2)
      SPM(IM1) = ANS
      USP(IM1)= BV(S1)
C
C  MOVE ALONG HORIZONTAL MESH LINE UNTIL MESH LINE INTERSECTS BOUNDARY
C
   80 ITI= IT
   90 ITV2 = ITV(MBO,3)
      IF (MBI2.GT.MBO) ITV2 = ITV(MBOP1,3)+NBBI
      IF (IM.EQ.MBOP1.AND.IT.GE.ITV2.AND.IT.LE.ITV(MBO,2)) ADD=.TRUE.
      IF (ADD) ADDL = .TRUE.
      IF (ADD) ITI=IT-NBBI
      IF (ADD.AND.S1.EQ.3) USP(IM1)=BV(S1)+1.
      IF (IM.LT.MBI.OR.IM.GT.MBO2) GO TO 120
      DO 100 SURF=1,3,2
      IF (IM.LE.MBI2.AND.SURF.EQ.3) GO TO 100
      IF (IM.EQ.IM2.AND.S1.EQ.4.AND.SURF.EQ.3) GO TO 100
      ITVIM1 = ITV(IM-1,SURF)
      IF (IM.EQ.MBOP1.AND.ADD) ITVIM1 = ITVIM1-NBBI
      IF (ITI.LT.ITV(IM,SURF).AND.ITI.GE.ITVIM1) GO TO 140
  100 CONTINUE
      SURF = 3
      IF(IM.EQ.MBI2.AND.IT.EQ.ITV(MBI2,3)-1.AND.ITV(MBI2,3)-ITV(MBI2,4)
     1.EQ.2) GO TO 140
      DO 110 SURF=2,4,2
      IF (IM.GT.MBO.AND.SURF.EQ.2) GO TO 110
      IF (IM.EQ.IM2.AND.S1.EQ.3.AND.SURF.EQ.4) GO TO 110
      ITVIM1 = ITV(IM-1,SURF)
      IF (IM.EQ.MBOP1.AND.ADD) ITVIM1 = ITVIM1-NBBI
      IF (ITI.GT.ITV(IM,SURF).AND.ITI.LE.ITVIM1) GO TO 140
  110 CONTINUE
  120 SPM(IM) = MV(IM)
      IP = IPF(IM,ITI)
      USP(IM) = U(IP)
      IF (ADD) USP(IM) = USP(IM)+1.
      IF (IM.EQ.MM) GO TO 130
      IM= IM+1
      GO TO 90
C
C  FINAL POINT IS ON BOUNDARY G-H
C
  130 IMT = MM
      GO TO 150
C
C  FINAL POINT IS ON A BLADE SURFACE
C
  140 ST = SURF
      IMT=IM
      IMTM1= IMT-1
      TH = FLOAT(ITI)*HT
      IF (ST.EQ.3.AND.IMT.LE.MBO) TH = TH-FLOAT(NBBI)*HT
      MVIM1 = MV(IMTM1)
      IF ((IMTM1.EQ.MBI).AND.(ST.EQ.1.OR.ST.EQ.2))
     1   MVIM1 = MVIM1+(MV(IMT)-MVIM1)/1000.
      IF ((IMTM1.EQ.MBI2).AND.(ST.EQ.3.OR.ST.EQ.4).AND.
     1   (ITV(MBI2,3)-ITV(MBI2,4).EQ.2.OR.ITV(MBI2,4)-ITV(MBI2,3).EQ.
     2   NBBI-2)) MVIM1 = MVIM1+(MV(IMT)-MVIM1)/1000.
      LER(2)=14
C     BLCD (VIA ROOT) CALL NO. 14
      IF(ST.EQ.1.AND.IMT.NE.MBI)CALL ROOT(MVIM1,MV(IMT),TH,BL1,
     1DTLR,ANS,AAA)
      LER(2)=15
C     BLCD (VIA ROOT) CALL NO. 15
      IF(ST.EQ.3.AND.IMT.NE.MBI2)CALL ROOT(MVIM1,MV(IMT),TH,BL3,
     1DTLR,ANS,AAA)
      LER(2)=16
C     BLCD (VIA ROOT) CALL NO. 16
      IF(ST.EQ.2)CALL ROOT(MVIM1,MV(IMT),TH,BL2,DTLR,ANS,AAA)
      LER(2)=17
C     BLCD (VIA ROOT) CALL NO. 17
      IF(ST.EQ.4)CALL ROOT(MVIM1,MV(IMT),TH,BL4,DTLR,ANS,AAA)
      IF(ST.EQ.1.AND.IMT.EQ.MBI) ANS = MV(MBI)
      IF(ST.EQ.3.AND.IMT.EQ.MBI2) ANS = MV(MBI2)
      SPM(IMT) = ANS
      USP(IMT)= BV(ST)
      IF(ST.EQ.3.AND.IMT.LE.MBO) USP(IMT) = BV(4)
      IF(ADD) USP(IMT) = USP(IMT)+1.
      IF (S1.EQ.3.AND.ST.EQ.2) USP(IM1) = BV(S1)+1.
      IF (S1.EQ.3.AND.ST.EQ.3) USP(IM1) = BV(S1)+1.
C
C  CALCULATE RHO*W-SUB-THETA AND THEN RHO*W AND BETA IN THE REGION
C
  150 NSP= IMT-IM1+1
      CALL SPLINE(SPM(IM1),USP(IM1),NSP,DUDM(IM1),AAA(IM1))
      FIRST=1
      IF (IM1.NE.1) FIRST=IM2
      LAST= MM
      IF (IMT.NE.MM) LAST=IMTM1
      IF (FIRST.GT.LAST) GO TO 170
      ITI = IT
      IF (FIRST.GT.MBOP1.AND.ADD) ITI=IT-NBBI
      DO 160 I=FIRST,LAST
      IF (ADD.AND.I.EQ.MBOP1) ITI=IT-NBBI
      RWT = -DUDM(I)*WTFL/BE(I)
      IP = IPF(I,ITI)
      W(IP) = SQRT(RWT**2+RWM(IP)**2)
  160 BETA(IP) = ATAN2(RWT,RWM(IP))*57.295779
C
C  CALCULATE RHO*W ON THE BLADE SURFACES
C
  170 IF (IM1.EQ.1) GO TO 180
      CALL SEARCH (SPM(IM1),S1,IHS)
      ANS = -DUDM(IM1)*WTFL/BEH(IHS,S1)
      WTB(IHS,S1) = ABS(ANS)*SQRT(1.+1./(RMH(IHS,S1)*DTDMH(IHS,S1))**2)
  180 IF(IMT.EQ.MM) GO TO 200
      CALL SEARCH (SPM(IMT),ST,IHS)
      ANS = -DUDM(IMT)*WTFL/BEH(IHS,ST)
      WTB(IHS,ST) = ABS(ANS)*SQRT(1.+1./(RMH(IHS,ST)*DTDMH(IHS,ST))**2)
  190 GO TO 20
  200 IT = IT+1
      GO TO 10
      END
$IBFTC SEARCH  DEBUG
      SUBROUTINE SEARCH (DIST,SURF,IS)
C
C  SEARCH LOCATES THE POSITION OF A GIVEN VALUE OF M IN THE MH ARRAY
C
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      DO 10 I=1,100
      IF (ABS(MH(I,SURF)-DIST).GT.DMLR) GO TO 10
      IS = I
      RETURN
   10 CONTINUE
      WRITE (6,1000) DIST,SURF
      STOP
 1000 FORMAT (38HL SEARCH CANNOT FIND M IN THE MH ARRAY/7H DIST =,G14.6,
     110X,6HSURF =,G14.6)
      END
$IBFTC VELOCY  DEBUG
      SUBROUTINE VELOCY
C
C  VELOCY CALLS SUBROUTINES TO CALCULATE DENSITIES AND VELOCITIES
C  THROUGHOUT THE REGION AND ON THE BLADE SURFACES, AND IT PLOTS
C  THE SURFACE VELOCITIES
C
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      DIMENSION KKK(18)
      DIMENSION W(2000),RWM(2000),BETA(2000),WMB(100,4),WTB(100,4),
     1XDOWN(800),YACROS(800)
      EQUIVALENCE (A(1,1),W(1)),(A(1,2),RWM(1)),(A(1,3),BETA(1)),
     1(A(1,4),WMB(1)),(A(401,4),WTB(1)),(A(801,4),XDOWN(1)),
     2(K(1),YACROS(1))
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      DATA KKK(4)/1H*/,KKK(6)/1H0/,KKK(8)/1H=/,KKK(10)/1H(/,
     1KKK(12)/1H+/,KKK(14)/1HX/,KKK(16)/1H$/,KKK(18)/1H)/
C
C  CALL VELP, VELBB, AND VELSUR THROUGHOUT THE REGION
C
      CALL VELP(1,MBIM1,1,2)
      IF (MBI2.GT.MBO) GO TO 10
      CALL VELBB(MBI,MBI2M1,1,2)
      CALL VELBB(MBI2,MBO,1,4)
      CALL VELBB(MBI2,MBO,3,2)
      CALL VELBB(MBOP1,MBO2,3,4)
      GO TO 20
   10 CALL VELBB(MBI,MBO,1,2)
      CALL VELP(MBOP1,MBI2M1,3,4)
      CALL VELBB(MBI2,MBO2,3,4)
   20 CALL VELP(MBO2P1,MM,3,4)
      CALL VELSUR
C
C  PREPARE INPUT ARRAYS FOR PLOT OF VELOCITIES
C
      IF (SURVL.LE.0) RETURN
      NP2 = 0
C   SURFACES 1 TO 4 - TANGENTIAL COMPONENTS
      DO 50 SURF=1,4
      NP1 = NP2
      IMSS = IMS(SURF)
      IF(IMSS.LT.1) GO TO 40
      DO 30 IHS=1,IMSS
      IF (ABS(DTDMH(IHS,SURF)*RMH(IHS,SURF)).LT..57735) GO TO 30
      NP1 = NP1+1
      YACROS(NP1) = WTB(IHS,SURF)
      XDOWN(NP1) = MH(IHS,SURF)
   30 CONTINUE
   40 KKK(2*SURF+1) = NP1-NP2
   50 NP2 = NP1
C   SURFACES 1 AND 2 - MERIDIONAL COMPONENTS
      DO 80 SURF=1,2
      NP1 = NP2
      IF (MBIP1.GT.MBOM1) GO TO 70
      DO 60 IM=MBIP1,MBOM1
      IF (ABS(DTDMV(IM,SURF)*RM(IM)).GT.1.7321) GO TO 60
      NP1 = NP1+1
      YACROS(NP1) = WMB(IM,SURF)
      XDOWN(NP1) = MV(IM)
   60 CONTINUE
   70 KKK(2*SURF+9) = NP1-NP2
   80 NP2 = NP1
C   SURFACES 3 AND 4 - MERIDIONAL COMPONENTS
      DO 110 SURF=3,4
      NP1 = NP2
      IF(MBI2P1.GT.MBO2M1) GO TO 100
      DO 90 IM=MBI2P1,MBO2M1
      IF (ABS(DTDMV(IM,SURF)*RM(IM)).GT.1.7321) GO TO 90
      NP1 = NP1+1
      YACROS(NP1) = WMB(IM,SURF)
      XDOWN(NP1) = MV(IM)
   90 CONTINUE
  100 KKK(2*SURF+9) = NP1-NP2
  110 NP2 = NP1
C
C  PLOT VELOCITIES
C
      KKK(1) = 1
      KKK(2) = 8
      P = 5.
      WRITE(6,1000)
      CALL PLOTMY(XDOWN,YACROS,KKK,P)
      WRITE(6,1010)
      RETURN
 1000 FORMAT(2HPT,50X,24HBLADE SURFACE VELOCITIES)
 1010 FORMAT (2HPL,37X,63HVELOCITY(W) VS. MERIDIONAL STREAMLINE DISTANCE
     1(M) DOWN THE PAGE  /2HPL/
     22HPL,50X,50H+ - BLADE SURFACE 1, BASED ON MERIDIONAL COMPONENT/
     32HPL,50X,50H* - BLADE SURFACE 1, BASED ON TANGENTIAL COMPONENT/
     42HPL,50X,50HX - BLADE SURFACE 2, BASED ON MERIDIONAL COMPONENT/
     52HPL,50X,50H0 - BLADE SURFACE 2, BASED ON TANGENTIAL COMPONENT/
     62HPL,50X,50H$ - BLADE SURFACE 3, BASED ON MERIDIONAL COMPONENT/
     72HPL,50X,50H= - BLADE SURFACE 3, BASED ON TANGENTIAL COMPONENT/
     82HPL,50X,50H) - BLADE SURFACE 4, BASED ON MERIDIONAL COMPONENT/
     92HPL,50X,50H( - BLADE SURFACE 4, BASED ON TANGENTIAL COMPONENT)
      END
$IBFTC VEL     DEBUG
      SUBROUTINE VEL
C
C  VEL CALCULATES DENSITIES AND VELOCITIES FROM THE PRODUCT OF
C  DENSITY TIMES VELOCITY
C
      COMMON SRW,ITER,IEND,LER(2),NER(2)
      COMMON /AUKRHO/ A(2000,4),U(2000),K(2000),RHO(2000)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /RHOS/ RHOHB(100,4),RHOVB(100,4)
      DIMENSION WWCRM(100,4),WWCRT(100,4),SURFL(100,4)
      DIMENSION W(2000),RWM(2000),BETA(2000),WMB(100,4),WTB(100,4),
     1XDOWN(800),YACROS(800)
      EQUIVALENCE (A(1,1),W(1)),(A(1,2),RWM(1)),(A(1,3),BETA(1)),
     1(A(1,4),WMB(1)),(A(401,4),WTB(1)),(A(801,4),XDOWN(1)),
     2(K(1),YACROS(1))
C
C  VELP CALCULATES ALONG VERTICAL MESH LINES WHICH DO NOT
C  INTERSECT BLADES
C
      ENTRY VELP(FIRST,LAST,UPPER,LOWER)
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      IF (FIRST.GT.LAST) RETURN
      IF (FIRST.EQ.1.AND.INTVL.GT.0) WRITE(6,1000)
      IF (FIRST.EQ.1) RELER = .0
      DO 20 IM=FIRST,LAST
      IPU = IV(IM)
      IPL = IPU+NBBI-1
      TWLMR = 2.*OMEGA*LAMBDA-(OMEGA*RM(IM))**2
      LER(1)=4
      DO 10 IP=IPU,IPL
C     DENSTY CALL NO. 4
      CALL DENSTY(W(IP),RHO(IP),ANS,TWLMR,CPTIP,EXPON,RHOIP,GAM,AR,TIP)
   10 W(IP) = ANS
      IF (INTVL.LE.0) GO TO 20
      WRITE (6,1010) IM,(W(IP),BETA(IP),IP=IPU,IPL)
   20 CONTINUE
      RETURN
C
C  VELBB CALCULATES ALONG VERTICAL MESH LINES WHICH INTERSECT BLADES
C
      ENTRY VELBB(FIRST,LAST,UPPER,LOWER)
      IF (FIRST.GT.LAST) RETURN
      IF (FIRST.NE.MBI) GO TO 30
      SURFL(MBI,1) = 0.
      SURFL(MBI,2) = 0.
      SURFL(MBI2,3) = 0.
      SURFL(MBI2,4) = 0.
   30 DO 70 IM=FIRST,LAST
      ITVU = ITV(IM,UPPER)
      ITVL = ITV(IM,LOWER)
      IPUP1 = IPF(IM,ITVU)
      IPLM1 = IPF(IM,ITVL)
      TWLMR = 2.*OMEGA*LAMBDA-(OMEGA*RM(IM))**2
      WCR = SQRT(TGROG*TIP*(1.-TWLMR/CPTIP))
      IF (ITVL.LT.ITVU) GO TO 50
C  ALONG THE LINE BETWEEN BLADES
      LER(1)=5
      DO 40 IP=IPUP1,IPLM1
C     DENSTY CALL NO. 5
      CALL DENSTY(W(IP),RHO(IP),ANS,TWLMR,CPTIP,EXPON,RHOIP,GAM,AR,TIP)
   40 W(IP) = ANS
      IF (INTVL.LE.0) GO TO 50
      WRITE (6,1010) IM,(W(IP),BETA(IP),IP=IPUP1,IPLM1)
C  ON THE UPPER SURFACE
   50 RHOB = RHOVB(IM,UPPER)
      LER(1)=6
C     DENSTY CALL NO. 6
      CALL DENSTY(WMB(IM,UPPER),RHOVB(IM,UPPER),ANS,TWLMR,CPTIP,
     1EXPON,RHOIP,GAM,AR,TIP)
      WMB(IM,UPPER) = ANS
      WWCRM(IM,UPPER) = WMB(IM,UPPER)/WCR
      IF (IM.EQ.MBI.OR.(IM.EQ.MBI2.AND.UPPER.EQ.3)) GO TO 60
      DELTV = TV(IM-1,UPPER)-TV(IM,UPPER)
      IF (IM.EQ.MBOP1.AND.UPPER.EQ.3) DELTV = DELTV-PITCH
      SURFL(IM,UPPER) = SURFL(IM-1,UPPER) + SQRT((MV(IM)-MV(IM-1))**2 +
     1(DELTV*(RM(IM)+RM(IM-1))/2.)**2)
   60 RELER = AMAX1(RELER,ABS((RHOB-RHOVB(IM,UPPER))/RHOVB(IM,UPPER)))
C  ON THE LOWER SURFACE
      RHOB = RHOVB(IM,LOWER)
      LER(1)=7
C     DENSTY CALL NO. 7
      CALL DENSTY(WMB(IM,LOWER),RHOVB(IM,LOWER),ANS,TWLMR,CPTIP,
     1EXPON,RHOIP,GAM,AR,TIP)
      WMB(IM,LOWER) = ANS
      WWCRM(IM,LOWER) = WMB(IM,LOWER)/WCR
      IF (IM.EQ.MBI.OR.(IM.EQ.MBI2.AND.LOWER.EQ.4)) GO TO 70
      DELTV = TV(IM-1,LOWER)-TV(IM,LOWER)
      SURFL(IM,LOWER) = SURFL(IM-1,LOWER) + SQRT((MV(IM)-MV(IM-1))**2 +
     1(DELTV*(RM(IM)+RM(IM-1))/2.)**2)
   70 RELER = AMAX1(RELER,ABS((RHOB-RHOVB(IM,LOWER))/RHOVB(IM,LOWER)))
      RETURN
C
C  VELSUR CALCULATES ALONG A BLADE SURFACE
C
      ENTRY VELSUR
      DO 90 SURF=1,4
      IMSS = IMS(SURF)
      IF(IMSS.EQ.0) GO TO 90
      DO 80 IHS=1,IMSS
      TWLMR = 2.*OMEGA*LAMBDA-(OMEGA*RMH(IHS,SURF))**2
      WCR = SQRT(TGROG*TIP*(1.-TWLMR/CPTIP))
      RHOB = RHOHB(IHS,SURF)
      LER(1)=8
C     DENSTY CALL NO. 8
      CALL DENSTY(WTB(IHS,SURF),RHOHB(IHS,SURF),ANS,TWLMR,CPTIP,
     1EXPON,RHOIP,GAM,AR,TIP)
      WTB(IHS,SURF) = ANS
      WWCRT(IHS,SURF) = WTB(IHS,SURF)/WCR
   80 RELER = AMAX1(RELER,ABS((RHOB-RHOHB(IHS,SURF))/RHOHB(IHS,SURF)))
   90 CONTINUE
      IF (RELER.LT..001) IEND = IEND+1
      WRITE(6,1080) ITER,RELER
C
C   WRITE ALL BLADE SURFACE VELOCITIES
C
      IF (SURVL.LE.0) RETURN
      WRITE(6,1020)
      WRITE(6,1040) (MV(IM),WMB(IM,1),BETAV(IM,1),SURFL(IM,1),
     1WWCRM(IM,1),WMB(IM,2),BETAV(IM,2),SURFL(IM,2),WWCRM(IM,2),
     2IM=MBI,MBO)
      WRITE(6,1030)
      WRITE(6,1040) (MV(IM),WMB(IM,3),BETAV(IM,3),SURFL(IM,3),
     1WWCRM(IM,3),WMB(IM,4),BETAV(IM,4),SURFL(IM,4),WWCRM(IM,4),
     2IM=MBI2,MBO2)
      WRITE(6,1050)
      DO 100 SURF=1,4
      IMSS = IMS(SURF)
      IF(IMSS.LT.1) GO TO 100
      WRITE(6,1060) SURF
      WRITE(6,1070) (MH(IHS,SURF),WTB(IHS,SURF),BETAH(IHS,SURF),
     1WWCRT(IHS,SURF),IHS=1,IMSS)
  100 CONTINUE
      RETURN
 1000 FORMAT(1H1////40X,34HVELOCITIES AT INTERIOR MESH POINTS//)
 1010 FORMAT(1HL,3HIM=,I3,5(24H   VELOCITY   ANGLE(DEG))/
     1(5X,5(G15.4,F9.2)))
 1020 FORMAT(1H1////16X,1H*,18X,49HSURFACE VELOCITIES BASED ON MERIDIONA
     1L COMPONENTS,43X,1H*/16X,1H*,53X,1H*,56X,1H*/16X,1H*,19X,15HBLADE
     2SURFACE 1,19X,1H*,20X,15HBLADE SURFACE 2,21X,1H*/7X,1HM,8X,1H*,2(3
     3X,8HVELOCITY,3X,23HANGLE(DEG) SURF. LENGTH,5X,5HW/WCR,6X,1H*,3X))
 1030 FORMAT(///16X,1H*,19X,15HBLADE SURFACE 3,19X,1H*,20X,15HBLADE SURF
     1ACE 4,21X,1H*/7X,1HM,8X,1H*,2(3X,8HVELOCITY,3X,23HANGLE(DEG) SURF.
     2 LENGTH,5X,5HW/WCR,6X,1H*,3X))
 1040 FORMAT(1H ,G13.4,3H  *,G12.4,F9.2,2G15.4,6H  *   ,G12.4,F9.2,
     12G15.4,3H  *)
 1050 FORMAT(1H1////3X,49HSURFACE VELOCITIES BASED ON TANGENTIAL COMPONE
     1NTS)
 1060 FORMAT(//22X,15HBLADE SURFACE  ,I1/7X,1HM,10X,8HVELOCITY,3X,10HANG
     1LE(DEG),3X,5HW/WCR)
 1070 FORMAT(1H ,2G13.4,F9.2,G15.4)
 1080 FORMAT(14HLITERATION NO.,I3,3X,36HMAXIMUM RELATIVE CHANGE IN DENSI
     1TY =,G11.4)
      END
$IBFTC SPLINE  DEBUG
      SUBROUTINE SPLINE (X,Y,N,SLOPE,EM)
C
C  SPLINE CALCULATES FIRST AND SECOND DERIVATIVES AT SPLINE POINTS
C  END CONDITION - SECOND DERIVATIVES ARE THE SAME AT END POINT AND
C  ADJACENT POINT
C
      COMMON Q/BOX/S(100),A(100),B(100),C(100),F(100),W(100),SB(100),
     1G(200)
      DIMENSION X(N),Y(N),EM(N),SLOPE(N)
      INTEGER Q
      DO 10 I=2,N
   10 S(I)=X(I)-X(I-1)
      NO=N-1
      IF(NO.LT.2) GO TO 30
      DO 20 I=2,NO
      A(I)=S(I)/6.
      B(I)=(S(I)+S(I+1))/3.
      C(I)=S(I+1)/6.
   20 F(I)=(Y(I+1)-Y(I))/S(I+1)-(Y(I)-Y(I-1))/S(I)
   30 A(N) = -.5
      B(1)=1.
      B(N)=1.
      C(1) = -.5
      F(1)=0.
      F(N)=0.
      W(1)=B(1)
      SB(1)=C(1)/W(1)
      G(1)=0.
      DO 40 I=2,N
      W(I)=B(I)-A(I)*SB(I-1)
      SB(I)=C(I)/W(I)
   40 G(I)=(F(I)-A(I)*G(I-1))/W(I)
      EM(N)=G(N)
      DO 50 I=2,N
      K=N+1-I
   50 EM(K)=G(K)-SB(K)*EM(K+1)
      SLOPE(1)=-S(2)/6.*(2.*EM(1)+EM(2))+(Y(2)-Y(1))/S(2)
      DO 60 I=2,N
   60 SLOPE(I)=S(I)/6.*(2.*EM(I)+EM(I-1))+(Y(I)-Y(I-1))/S(I)
      IF (Q.EQ.13) WRITE(6,1000) N,(X(I),Y(I),SLOPE(I),EM(I),I=1,N)
      RETURN
 1000 FORMAT (2X,15HNO. OF POINTS =,I3/10X,1HX,19X,1HY,19X,5HSLOPE,15X,
     12HEM/(4F20.8))
      END
$IBFTC SPLINT  DEBUG
      SUBROUTINE SPLINT (X,Y,N,Z,MAX,YINT,DYDX)
C
C  SPLINT CALCULATES INTERPOLATED POINTS AND DERIVATIVES
C  FOR A SPLINE CURVE
C  END CONDITION - SECOND DERIVATIVES ARE THE SAME AT END POINT AND
C  ADJACENT POINT
C
      COMMON Q/BOX/S(100),A(100),B(100),C(100),F(100),W(100),SB(100),
     1G(100),EM(100)
      DIMENSION   X(N),Y(N),Z(MAX),YINT(MAX),DYDX(MAX)
      INTEGER Q
      IF(MAX.LE.0) RETURN
      III = Q
      DO 10 I=2,N
   10 S(I)=X(I)-X(I-1)
      NO=N-1
      IF(NO.LT.2) GO TO 30
      DO 20 I=2,NO
      A(I)=S(I)/6.0
      B(I)=(S(I)+S(I+1))/3.0
      C(I)=S(I+1)/6.0
   20 F(I)=(Y(I+1)-Y(I))/S(I+1)-(Y(I)-Y(I-1))/S(I)
   30 A(N) = -.5
      B(1)=1.0
      B(N)=1.0
      C(1) = -.5
      F(1)=0.0
      F(N)=0.0
      W(1)=B(1)
      SB(1)=C(1)/W(1)
      G(1)=0.0
      DO 40 I=2,N
      W(I)=B(I)-A(I)*SB(I-1)
      SB(I)=C(I)/W(I)
   40 G(I)=(F(I)-A(I)*G(I-1))/W(I)
      EM(N)=G(N)
      DO 50 I=2,N
      K=N+1-I
   50 EM(K)=G(K)-SB(K)*EM(K+1)
      DO 140 I=1,MAX
      K=2
      IF(Z(I)-X(1)) 70,60,90
   60 YINT(I)=Y(1)
      GO TO 130
   70 IF(Z(I).GE.(1.1*X(1)-.1*X(2))) GO TO 120
      WRITE (6,1000) Z(I)
      Q = 16
      GO TO 120
   80 K=N
      IF(Z(I).LE.(1.1*X(N)-.1*X(N-1))) GO TO 120
      WRITE (6,1000) Z(I)
      Q = 16
      GO TO 120
   90 IF(Z(I)-X(K)) 120,100,110
  100 YINT(I)=Y(K)
      GO TO 130
  110 K=K+1
      IF(K-N) 90,90,80
  120 YINT(I) = EM(K-1)*(X(K)-Z(I))**3/6./S(K)+EM(K)*(Z(I)-X(K-1))**3/6.
     1/S(K)+(Y(K)/S(K)-EM(K)*S(K)/6.)*(Z(I)-X(K-1))+(Y(K-1)/S(K)-EM(K-1)
     2*S(K)/6.)*(X(K)-Z(I))
  130 DYDX(I)=-EM(K-1)*(X(K)-Z(I))**2/2.0/S(K)+EM(K)*(X(K-1)-Z(I))**2/2.
     10/S(K)+(Y(K)-Y(K-1))/S(K)-(EM(K)-EM(K-1))*S(K)/6.0
  140 CONTINUE
      MXA = MAX0(N,MAX)
      IF(Q.EQ.16) WRITE(6,1010) N,MAX,(X(I),Y(I),Z(I),YINT(I),DYDX(I),
     1I=1,MXA)
      Q = III
      RETURN
 1000 FORMAT (54H SPLINT USED FOR EXTRAPOLATION.  EXTRAPOLATED VALUE = ,
     1G14.6)
 1010 FORMAT (2X,21HNO. OF POINTS GIVEN =,I3,30H, NO. OF INTERPOLATED PO
     1INTS =,I3/10X,1HX,19X,1HY,16X,11HX-INTERPOL.,9X,11HY-INTERPOL.,
     28X,14HDYDX-INTERPOL./(5E20.8))
      END
$IBFTC SPLN22  DEBUG
      SUBROUTINE SPLN22 (X,Y,Y1P,YNP,N,SLOPE,EM)
C
C  SPLN22 CALCULATES FIRST AND SECOND DERIVATIVES AT SPLINE POINTS
C  END CONDITION - DERIVATIVES SPECIFIED AT END POINTS
C
      COMMON Q/BOX/S(100),A(100),B(100),C(100),F(100),W(100),SB(100),
     1G(200)
      DIMENSION X(N),Y(N),EM(N),SLOPE(N)
      INTEGER Q
      DO 10 I=2,N
   10 S(I)=X(I)-X(I-1)
      NO=N-1
      IF(NO.LT.2) GO TO 30
      DO 20 I=2,NO
      A(I)=S(I)/6.
      B(I)=(S(I)+S(I+1))/3.
      C(I)=S(I+1)/6.
   20 F(I)=(Y(I+1)-Y(I))/S(I+1)-(Y(I)-Y(I-1))/S(I)
   30 A(N) = S(N)/6.
      B(1)=S(2)/3.
      B(N) = S(N)/3.
      C(1)=S(2)/6.
      F(1)=(Y(2)-Y(1))/S(2)-Y1P
      F(N) = YNP-(Y(N)-Y(N-1))/S(N)
      W(1)=B(1)
      SB(1)=C(1)/W(1)
      G(1)=F(1)/W(1)
      DO 40 I=2,N
      W(I)=B(I)-A(I)*SB(I-1)
      SB(I)=C(I)/W(I)
   40 G(I)=(F(I)-A(I)*G(I-1))/W(I)
      EM(N)=G(N)
      DO 50 I=2,N
      K=N+1-I
   50 EM(K)=G(K)-SB(K)*EM(K+1)
      SLOPE(1)=-S(2)/6.*(2.*EM(1)+EM(2))+(Y(2)-Y(1))/S(2)
      DO 60 I=2,N
   60 SLOPE(I)=S(I)/6.*(2.*EM(I)+EM(I-1))+(Y(I)-Y(I-1))/S(I)
      IF (Q.EQ.18) WRITE(6,1000) N,(X(I),Y(I),SLOPE(I),EM(I),I=1,N)
      RETURN
 1000 FORMAT (2X,15HNO. OF POINTS =,I3/10X,1HX,19X,1HY,19X,5HSLOPE,15X,
     12HEM/(4F20.8))
      END
$IBFTC ROOT    DEBUG
      SUBROUTINE ROOT(A,B,Y,FUNCT,TOLERY,X,DFX)
C
C  ROOT FINDS A ROOT FOR (FUNCT MINUS Y) IN THE INTERVAL (A,B)
C
      COMMON SRW,ITER,IEND,LER(2),NER(2)
      INTEGER SRW
      IF (SRW.EQ.21) WRITE(6,1000) A,B,Y,TOLERY
      TOLERX= (B-A)/1000.
      AB2= (A+B)/2.
      I = 0
      X = A
   10 CALL FUNCT(X,FX,DFX,INF)
      IF (SRW.EQ.21) WRITE(6,1010) I,X,FX,DFX,INF
      IF (ABS(Y-FX).LT.TOLERY) RETURN
      IF (I.GE.1000) GO TO 30
      I= I+1
      IF (INF.NE.0 .OR. DFX.EQ.0.) GO TO 20
      X= (Y-FX)/DFX+X
      IF (X.GE.A .AND. X.LE.B) GO TO 10
      X = A+TOLERX*FLOAT(I)
      IF(I.EQ.1) X = B
      GO TO 10
   20 IF (X.LT.AB2) X=X+TOLERX
      IF (X.GE.AB2) X=X-TOLERX
      GO TO 10
   30 WRITE(6,1020) LER(2),A,B,Y
      STOP
 1000 FORMAT (32H1INPUT ARGUMENTS FOR ROOT -- A =G13.5,3X,3HB =,G13.5,
     13X,3HY =,G13.5,3X,8HTOLERY =,G13.5/17H  ITER. NO.     X,17X,
     22HFX,15X,3HDFX,10X,3HINF)
 1010 FORMAT (5X,I3,G16.5,2G18.5,I6)
 1020 FORMAT (14HLROOT CALL NO.,I3/47H ROOT HAS FAILED TO CONVERGE IN 10
     100 ITERATIONS/4H A =,G14.6,10X,3HB =,G14.6,10X,3HY =,G14.6)
      END
$IBFTC BLCD    DEBUG
      SUBROUTINE BLCD
C
C  BLCD CALCULATES BLADE THETA COORDINATE AS A FUNCTION OF M
C
      COMMON SRW,ITER,IEND,LER(2),NER(2)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      COMMON /GEOMIN/ CHORD(4),STGR(4),MLE(4),THLE(4),RMI(4),RMO(4),
     1RI(4),RO(4),BETI(4),BETO(4),NSPI(4),MSP(50,4),THSP(50,4)
      COMMON /BLCDCM/ EM(50,4),INIT(4)
      ENTRY  BL1(M,THETA,DTDM,INF)
      INTEGER BLDAT,AANDK,ERSOR,STRFN,SLCRD,SURVL,AATEMP,SURF,SURFBV,
     1FIRST,UPPER,UPPRBV,S1,ST,SRW
      REAL K,KAK,LAMBDA,LMAX,MH,MLE,MR,MSL,MSP,MV,MVIM1
      REAL M,MMLE,MSPMM,MMMSP
      SURF= 1
      SIGN= 1.
      GO TO 10
      ENTRY  BL2(M,THETA,DTDM,INF)
      SURF= 2
      SIGN= -1.
      GO TO 10
      ENTRY  BL3(M,THETA,DTDM,INF)
      SURF= 3
      SIGN= 1.
      GO TO 10
      ENTRY  BL4(M,THETA,DTDM,INF)
      SURF= 4
      SIGN= -1.
   10 INF= 0
      NSP= NSPI(SURF)
      IF (INIT(SURF).EQ.13) GO TO 30
      INIT(SURF)= 13
C
C  INITIAL CALCULATION OF FIRST AND LAST SPLINE POINTS ON BLADE
C
      AA = BETI(SURF)/57.295779
      AA = SIN(AA)
      MSP(1,SURF) =  RI(SURF)*(1.-SIGN*AA)
      BB = SQRT(1.-AA**2)
      THSP(1,SURF) = SIGN*BB*RI(SURF)/RMI(SURF)
      BETI(SURF) = AA/BB/RMI(SURF)
      AA = BETO(SURF)/57.295779
      AA = SIN(AA)
      MSP(NSP,SURF) = CHORD(SURF)-RO(SURF)*(1.+SIGN*AA)
      BB = SQRT(1.-AA**2)
      THSP(NSP,SURF) = STGR(SURF)+SIGN*BB*RO(SURF)/RMO(SURF)
      BETO(SURF) = AA/BB/RMO(SURF)
      DO 20 IA=1,NSP
      MSP(IA,SURF)= MSP(IA,SURF)+MLE(SURF)
   20 THSP(IA,SURF)= THSP(IA,SURF)+THLE(SURF)
      CALL SPLN22(MSP(1,SURF),THSP(1,SURF),BETI(SURF),BETO(SURF),NSP,
     1 AAA,EM(1,SURF))
      IF(BLDAT.LE.0) GO TO 30
      IF (SURF.EQ.1) WRITE(6,1000)
      WRITE(6,1010) SURF
      WRITE (6,1020) (MSP(IA,SURF),THSP(IA,SURF),AAA(IA),EM(IA,SURF),
     1IA=1,NSP)
C
C  BLADE COORDINATE CALCULATION
C
   30 KK = 2
      IF (M.GT.MSP(1,SURF)) GO TO 50
C
C  AT LEADING EDGE RADIUS
C
      MMLE= M-MLE(SURF)
      IF (MMLE.LT.-DMLR) GO TO 90
      MMLE= AMAX1(0.,MMLE)
      THETA= SQRT(MMLE*(2.*RI(SURF)-MMLE))*SIGN
      IF (THETA.EQ.0.) GO TO 40
      RMM= RI(SURF)-MMLE
      DTDM= RMM/THETA/RMI(SURF)
      THETA= THETA/RMI(SURF)+THLE(SURF)
      RETURN
   40 INF= 1
      DTDM = 1.E10*SIGN
      THETA= THLE(SURF)
      RETURN
C
C  ALONG SPLINE CURVE
C
   50 IF (M.LE.MSP(KK,SURF)) GO TO 60
      IF (KK.GE.NSP) GO TO 70
      KK = KK+1
      GO TO 50
   60 S= MSP(KK,SURF)-MSP(KK-1,SURF)
      EMKM1= EM(KK-1,SURF)
      EMK= EM(KK,SURF)
      MSPMM= MSP(KK,SURF)-M
      MMMSP= M-MSP(KK-1,SURF)
      THK= THSP(KK,SURF)/S
      THKM1= THSP(KK-1,SURF)/S
      THETA= EMKM1*MSPMM**3/6./S + EMK*MMMSP**3/6./S + (THK-EMK*S/6.)*
     1MMMSP+(THKM1-EMKM1*S/6.)*MSPMM
      DTDM= -EMKM1*MSPMM**2/2./S + EMK*MMMSP**2/2./S + THK-THKM1-(EMK-
     1EMKM1)*S/6.
      RETURN
C
C  AT TRAILING EDGE RADIUS
C
   70 CMM= CHORD(SURF)+MLE(SURF)-M
      IF (CMM.LT.-DMLR) GO TO 90
      CMM= AMAX1(0.,CMM)
      THETA= SQRT(CMM*(2.*RO(SURF)-CMM))*SIGN
      IF (THETA.EQ.0.) GO TO 80
      RMM= RO(SURF)-CMM
      DTDM = -RMM/THETA/RMO(SURF)
      THETA = STGR(SURF)+THETA/RMO(SURF)+THLE(SURF)
      RETURN
   80 INF= 1
      DTDM = -1.E10*SIGN
      THETA= THLE(SURF)+STGR(SURF)
      RETURN
C
C  ERROR RETURN
C
   90 WRITE(6,1030) LER(2),M,SURF
      STOP
 1000 FORMAT (1H1,13X,33HBLADE DATA AT INPUT SPLINE POINTS)
 1010 FORMAT(1HL,17X,16HBLADE    SURFACE,I4)
 1020 FORMAT (7X ,1HM,10X,5HTHETA,10X,10HDERIVATIVE,5X,10H2ND DERIV. /
     1(4G15.5))
 1030 FORMAT (14HLBLCD CALL NO.,I3/33H M COORDINATE IS NOT WITHIN BLADE/
     14H M =,G14.6,10X,6HSURF =,G14.6)
      END
$IBFTC DENSTY  DEBUG
      SUBROUTINE DENSTY(RHOW,RHO,VEL,TWLMR,CPTIP,EXPON,RHOIP,GAM,AR,TIP)
C
C  DENSTY CALCULATES DENSITY AND VELOCITY FROM THE WEIGHT FLOW PARAMETER
C  DENSITY TIMES VELOCITY
C
      COMMON SRW,ITER,IEND,LER(2),NER(2)
      VEL = RHOW/RHO
      IF (VEL.NE.0.) GO TO 10
      RHO = RHOIP
      RETURN
   10 TTIP = 1.-(VEL**2+TWLMR)/CPTIP
      IF(TTIP.LT.0.) GO TO 30
      TEMP = TTIP**(EXPON-1.)
      RHOT = RHOIP*TEMP*TTIP
      RHOWP = -VEL**2/GAM*RHOIP/AR*TEMP/TIP+RHOT
      IF(RHOWP.LE.0.) GO TO 30
      VELNEW = VEL+(RHOW-RHOT*VEL)/RHOWP
      IF(ABS(VELNEW-VEL)/VELNEW.LT..0001) GO TO 20
      VEL = VELNEW
      GO TO 10
   20 VEL = VELNEW
      RHO = RHOW/VEL
      RETURN
   30 TGROG = 2.*GAM*AR/(GAM+1.)
      VEL = SQRT(TGROG*TIP*(1.-TWLMR/CPTIP))
      RHO = RHOIP*(1.-(VEL**2+TWLMR)/CPTIP)**EXPON
      RWMORW = RHOW/RHO/VEL
      NER(1) = NER(1)+1
      WRITE(6,1000) LER(1),NER(1),RWMORW
      IF(NER(1).EQ.50) STOP
      RETURN
 1000 FORMAT(16HLDENSTY CALL NO.,I3/9H NER(1) =,I3/10H RHO*W IS ,F7.4,
     134H TIMES THE MAXIMUM VALUE FOR RHO*W)
      END
$IBFTC IPF     DEBUG
      FUNCTION IPF(IM,IT)
      COMMON /INP/ GAM,AR,TIP,RHOIP,WTFL,WTFLSP,OMEGA,ORF,BETAI,BETAO,
     1MBI,MBO,MBI2,MBO2,MM,NBBI,NBL,NRSP,MR(50),RMSP(50),BESP(50),
     2BLDAT,AANDK,ERSOR,STRFN,SLCRD,INTVL,SURVL
      COMMON /CALCON/ MBIM1,MBIP1,MBOM1,MBOP1,MBI2M1,MBI2P1,MBO2M1,
     1MBO2P1,MMM1,HM1,HM2,HM3,HT,DTLR,DMLR,PITCH,CP,EXPON,TWW,CPTIP,
     2TGROG,TBI,TBO,LAMBDA,TWL,ITMIN,ITMAX,NIP,IMS(4),BV(4),MV(100),
     3IV(101),ITV(100,4),TV(100,4),DTDMV(100,4),BETAV(100,4),
     4MH(100,4),DTDMH(100,4),BETAH(100,4),RMH(100,4),BEH(100,4),
     5RM(100),BE(100),DBDM(100),SAL(100),AAA(100)
      IPF = IV(IM)+IT-ITV(IM,1)-ITV(IM,3)-10000
      IF(IM.LT.MBI2.OR.IM.GT.MBO) RETURN
      IPF = IV(IM)+IT-ITV(IM,1)
      IF(IT.GE.ITV(IM,3)) IPF = IPF-ITV(IM,3)+ITV(IM,4)+1
      RETURN
      END
$IBMAP TIME1
HEY    OPSYN   BTTE                                                            1
       ENTRY   SYSONE
       ENTRY   TIME1                                                           2
TIME1  SAVE    4                                                               3
       CLA   DEC12                                                             4
       ZET     READY         SEE IF 40 IS BUSY                                 5
       TRA     *-1             YES                                             6
       STO     READY             NO                                            7
       HEY                                                                     8
       ZET     READY         SEE IF 40 HAS FINISHED                            9
       TRA     *-1           NYET                                             10
       CLA     READY+1       OK                                               11
       ANA     TIME2                                                          12
       ORA     TIME3                                                          13
       FAD     TIME3                                                          14
       STO*    3,4                                                            15
       RETURN  TIME1                                                          16
TIME2  OCT     77777777                                                       17
TIME3  OCT     233000000000                                                   18
DEC12  PZE     1,0,19        40 CODE TO READ THE CLOCK                        19
READY  EQU     18                                                             20
SYSONE PZE     1
       END                                                                    21
$IBFTC PISTUG  ALTIO
       SUBROUTINE PISTUG(ARRAY)                                                1
       DIMENSION ARRAY(1)                                                      2
      COMMON/JOLO/N,F,DX,XYX,FORY,STUG,LABOUT,TONLY,KSW64,KPWR,KFD,TLINX       3
       LOGICAL XYX,FORY, STUG,TONLY                                            4
  126 X1 = ARRAY(1)                                                            5
  128 IF(XYX) GO TO 133                                                        6
  130 DO 132 J = 2,N                                                           7
  132 X1 = AMIN1(X1,ARRAY(J))                                                  8
  133 IF(STUG) X1=F                                                            9
  134 XN = 0.0                                                                10
  136 DO 146 J = 1,N                                                          11
  138 DIF = ABS(X1-ARRAY(J))                                                  12
  140 IF(DIF.LE.XN)GO TO 146                                                  13
  142 XN=DIF                                                                  14
  144 IHOLD=J                                                                 15
  146 CONTINUE                                                                16
  147 XN = ARRAY(IHOLD)                                                       17
  148 IF(KSW64.EQ.2) KPWR = KHAR(AMAX1(ABS(X1),ABS(XN)))                      18
  149 IF(TONLY) GO TO 240                                                     19
  150 TLIN=101.                                                               20
  152 IF(.NOT.FORY) TLIN = TLINX                                              21
  154 C5 = (XN-X1)/TLIN                                                       22
  156 C6 = ABS(C5)                                                            23
  158 IF(C6.EQ.0.) GO TO 300                                                  24
  159 K7 = KHAR(C6)                                                           25
  160 C8 = 10.**K7                                                            26
  162 C9 = C6/C8                                                              27
  164 IF((2.5-C9).LE.0.0) GO TO 172                                           28
  166 D=2.                                                                    29
  168 IF((2.0-C9).LE.0.) D=2.5                                                30
  170 GO TO 176                                                               31
  172 D=5.                                                                    32
  174 IF((5.-C9).LE.0.0) D=10.                                                33
  176 C11 = D*C8                                                              34
  178 DX = SIGN(C11,C5)                                                       35
  179 HUND = 100.*DX                                                          36
  240 K7 = KHAR(ABS(DX))                                                      37
  250 KFD = 0                                                                 38
  252 IF(K7) 260,270,254                                                      39
  254 IF(K7.GE.5) LABOUT=2                                                    40
  256 GO TO 270                                                               41
  260 KFD = 6                                                                 42
  262 IF(K7.LT.(-7)) LABOUT = 2                                               43
  264 IF(K7.GT.(-6)) KFD = -K7                                                44
  270 IF(STUG) GO TO 230                                                      45
  182 KC12 = INT(ABS(X1)/C11)                                                 46
  184 JJ = 1                                                                  47
  186 IF(X1) 188,192,190                                                      48
  188 JJ = 3                                                                  49
  190 IF(DX.LT.0.) JJ = JJ+1                                                  50
  192 GO TO (193,194,195,196),JJ                                              51
  193 KC14 = KC12                                                             52
       GO TO 198                                                              53
  194 KC14 = KC12+1                                                           54
       GO TO 198                                                              55
  195 KC14 = -KC12-1                                                          56
       GO TO 198                                                              57
  196 KC14 = -KC12                                                            58
  198 KC13 = MOD(KC12,10)                                                     59
       KC15 = KC12-KC13                                                       60
  199 KC18 = KC15                                                             61
  200 GO TO (212,202,202,210),JJ                                              62
  202 KC18 = KC18+10                                                          63
  204 IF(KC13.NE.9) GO TO 210                                                 64
  206 KC18 = KC14                                                             65
  208 GO TO 212                                                               66
  210 IF(X1.LT.0.)KC18 = -KC18                                                67
  212 F=C11*FLOAT(KC18)                                                       68
  214 IF(.NOT.FORY) GO TO 230                                                 69
  220 TEMP = F+HUND                                                           70
  222 GO TO (224,228,224,228),JJ                                              71
  224 IF(TEMP.GE.XN) GO TO 230                                                72
  226 GO TO 229                                                               73
  228 IF(TEMP.LE.XN) GO TO 230                                                74
  229 F=C11*FLOAT(KC14)                                                       75
  230 CONTINUE                                                                76
       RETURN                                                                 77
  300 DX=0.                                                                   78
       LABOUT=3                                                               79
       GO TO 230                                                              80
       END                                                                    81
$IBFTC PLOTMY  ALTIO
      SUBROUTINE PLOTMY(X,Y,K,P)                                               1
C                                                                              2
      COMMON/JOLO/N,F,DX,XYX,FORY,STUG,LABOUT,TONLY,KSW64,KPWR,KFD,TLINX       3
       LOGICAL XYX,FORY, STUG,TONLY,XGL,LS                                     4
      DIMENSION X(1),Y(1),P(1),K(1)                                            5
      DIMENSION FLS(3),FLAB(4),FYLAB(4),YLABEL(11),A(104),ELS(3)               6
       DIMENSION KPCSTD(6)                                                     7
      EQUIVALENCE (FLAB(3),IFLAB3),(FYLAB(3),IYLAB)                            8
       EQUIVALENCE (KPC,TPC)                                                   9
      DATA MASK1, MASK2,MASK4,MASK8,MASK16,MASK32,MASK64   /                  10
     1        O1,    O2,   O4,  O10,   O20,   O40,  O100   /                  11
      DATA FYLAB(1),FYLAB(2),FYLAB3,FYLAB(4)   /                              12
     1    6H(2HP ,,6H20X,11  ,6H F10.0,1H)  /                                 13
      DATA BLANK,XGRID,YGRID /1H ,1H-,1H1 /                                   14
       DATA RMARK/O726060606060/                                              15
       DATA(KPCSTD(I),I=1,6)/1H*,1H+,1H0,1HX,1H=,1HO/                         16
      DATA  FLS(1),FLS28,FLS264,FLS2B,FLS38,FLS3B /                           17
     1  6H(2HP+,,4H11X,,4H11X,,4HA11,,3HA6),6H F6.3)  /                       18
      DATA FLAB(1),FLAB(2),FLAB3,FLAB(4) /                                    19
     1   6H(2HP+,,4H17X,,6H  F9.0,1H)    /                                    20
C                                                                             21
  100 WRITE (6,500)                                                           22
  102 KODE=K(1)                                                               23
       KN=K(2)                                                                24
       NPTS=K(3)                                                              25
       LABOUT=1                                                               26
       FLAB(3) = FLAB3                                                        27
       LS = .FALSE.                                                           28
      FYLAB(3)=FYLAB3                                                         29
      KSW8=0                                                                  30
      KSW64=0                                                                 31
       KTL=1                                                                  32
       KSWI=1                                                                 33
       KSWII=1                                                                34
11000 IF(P(1)-2.5) 11002,11002,11010                                          35
11002 KSWI=2                                                                  36
11004 KTL=KN                                                                  37
11006 GO TO 11034                                                             38
11010 IF(P(1)-4.0)11020,11020,11012                                           39
11012 NPTST=0                                                                 40
11014 DO 11016 I=1,KN                                                         41
11016 NPTST=NPTST+K(2*I+1)                                                    42
11017 TLINX=55*(1+NPTST/(35*KN))                                              43
11018 GO TO112                                                                44
11020 KSWII=2                                                                 45
11021 KTIMES=KN-1                                                             46
11022 DO 11032 I=1,KTIMES                                                     47
11024 MM=I*NPTS                                                               48
11026 K(2*I+3)=NPTS                                                           49
11028 DO 11032 II=1,NPTS                                                      50
11030 L=MM+II                                                                 51
11032 Y(L)=Y(II)                                                              52
11034 NPTST=KN*NPTS                                                           53
11036 TLINX=55*(1+NPTS/35)                                                    54
  112 IF((AND(KODE,MASK1)).GT.0.)GO TO 114                                    55
11112 DO 11113 I=1,KN                                                         56
11113 K(2*I+2)=KPCSTD(I)                                                      57
  114 M=10                                                                    58
  116 IF((AND(KODE,MASK2)).GT.0.)M=P(3)                                       59
  117 IF (M.EQ.0)M=1000                                                       60
  118 NY=10                                                                   61
  120 IF((AND(KODE,MASK4)).GT.0.)NY = P(4)                                    62
  121 IF(NY.EQ.0) NY=1000                                                     63
  122 IF((AND(KODE,MASK64)).GT.0.)KSW64=2                                     64
  124 IF((AND(KODE,MASK8)).GT.0.)KSW8=1                                       65
C                                                                             66
  125 K864=KSW8+KSW64                                                         67
  126 IF(K864-2) 132,128,138                                                  68
  128 FLS(2)=FLS264                                                           69
  130 GO TO 139                                                               70
  132 FLS(2)=FLS28                                                            71
  134 FLS(3)=FLS38                                                            72
  136 GO TO 140                                                               73
  138 FLS(2)=FLS2B                                                            74
  139 FLS(3)=FLS3B                                                            75
C                                                                             76
  140 XYX=.FALSE.                                                             77
  142 FORY=.TRUE.                                                             78
  144  STUG=.FALSE.                                                           79
  146 TONLY=.FALSE.                                                           80
C                                                                             81
  148 IF((AND(KODE,MASK32).LE.0.))  GO TO 172                                 82
  151 STUG=.TRUE.                                                             83
  152 KSY=P(9)                                                                84
  154 PWR10Y=10.** (KSY-6)                                                    85
  156 FY =P(10)*PWR10Y                                                        86
  158 F = FY                                                                  87
C                                                                             88
  160 IF(P(5).GE.2.) GO TO 172                                                89
  162 TONLY=.TRUE.                                                            90
  164 DY= P(11)*PWR10Y                                                        91
  166 DX= DY                                                                  92
C                                                                             93
  172 N=NPTST                                                                 94
11172 CALL PISTUG(Y)                                                          95
  173 IF(DX.EQ.0.) GO TO 700                                                  96
  174 FY=F                                                                    97
  176 DY=DX                                                                   98
  190 IYLAB=IYLAB-KFD                                                         99
C                                                                            100
  200 XYX=.FALSE.                                                            101
  202 FORY=.FALSE.                                                           102
  204 STUG=.FALSE.                                                           103
  206 TONLY=.FALSE.                                                          104
C                                                                            105
  210 IF((AND(KODE,MASK16).LE.0.)) GO TO 232                                 106
  213 STUG=.TRUE.                                                            107
  214 KSX = P(6)                                                             108
  216 PWR10X=10.**(KSX-6)                                                    109
  218 FX= P(7)*PWR10X                                                        110
  220 F=FX                                                                   111
C                                                                            112
  222 IF(MOD(IFIX(P(5)),2).EQ.1) GO TO 232                                   113
  224 TONLY=.TRUE.                                                           114
  226 DX =P(8)*PWR10X                                                        115
  232 IF(KSWI.EQ.2)N=NPTS                                                    116
11233 ILIM=N                                                                 117
  233 CALL PISTUG(X)                                                         118
      IF(DX.EQ.0.) GO TO 700                                                 119
  234 FX=F                                                                   120
  248 IFLAB3=IFLAB3-KFD                                                      121
C                                                                            122
  250 IF(KSW64.EQ.0)GO TO 264                                                123
  252 KOUTX=-KPWR                                                            124
  256 F10X=10.**KPWR                                                         125
  260 WRITE (6,502) KOUTX                                                    126
C                                                                            127
  264 DO 278 I=1,11                                                          128
  266 TEMP = FY+FLOAT(I-1)*DY*10.                                            129
  268 ATEMP= ABS(TEMP)                                                       130
  270 IF (ATEMP.LT.1.E-7) TEMP = 0.                                          131
  272 IF (ATEMP.GE.1.E+7)LABOUT=2                                            132
  278 YLABEL(I)=TEMP                                                         133
  300 KSYLAB =1                                                              134
  302 WRITE (6,FYLAB) (YLABEL(I),I=1,11)                                     135
  304 GO TO (306,700),KSYLAB                                                 136
C                                                                            137
  306 KSYLAB =2                                                              138
  310 LCTR=0                                                                 139
      NCTR=1                                                                 140
      KOUT=1                                                                 141
      KQUIT= 1                                                               142
C                                                                            143
  320 IF(MOD(LCTR,M))328,322,328                                             144
  322 AFILL= XGRID                                                           145
  324 XGL =.TRUE.                                                            146
      GO TO 330                                                              147
  328 XGL =.FALSE.                                                           148
      AFILL=BLANK                                                            149
  330 DO 332 I=2,104                                                         150
  332 A(I) = AFILL                                                           151
  334 DO 336  I=2,104,NY                                                     152
  336 A(I)= YGRID                                                            153
      A(1) =BLANK                                                            154
  338 GO TO (340,400),KOUT                                                   155
C                                                                            156
  340 XMIN=1.E15                                                             157
11340 IMIN=1                                                                 158
11342 DO 11350 I=1,ILIM                                                      159
11344 IF(XMIN-X(I))11350,11350,11346                                         160
11346 XMIN=X(I)                                                              161
11348 IMIN=I                                                                 162
11350 CONTINUE                                                               163
  341 KX =(X(IMIN)-FX)/DX +.5                                                164
  342 IF(KX-LCTR)630,350,600                                                 165
  350 LS=.TRUE.                                                              166
11352 X(IMIN)=1.E15                                                          167
  372 J=1                                                                    168
  374 IF(KSW8.EQ.0) GO TO 380                                                169
  376  ELS(J)=P(IMIN+11)                                                     170
  378 J= J+1                                                                 171
  380 IF(KSW64.EQ.0) GO TO 11400                                             172
00382 ELS(J)=XMIN/F10X                                                       173
11400 DO 370 IM=1,KTL                                                        174
11402 LL=IMIN+(IM  -1)*NPTS                                                  175
11404 KY=(Y(LL)-FY)/DY+.5                                                    176
11420 IF(KSWI.EQ.2) GO TO 11440                                              177
11422 IK=0                                                                   178
       KLAST=2*KN+1                                                          179
11424 DO 11430 IL=3,KLAST,2                                                  180
11426 IK=IK+K(IL)                                                            181
11428 IF(IK-IMIN) 11430,11436,11436                                          182
11430 CONTINUE                                                               183
11432 LABOUT=5                                                               184
11434 GO TO 700                                                              185
11436 KPC=K(IL+1)                                                            186
11438 GO TO 353                                                              187
11440 KPC=K(2*IM+2)                                                          188
  353 KYL = KY+2                                                             189
  354 IF(KY.LT.0) GO TO 360                                                  190
  356 IF(KY.GT.101)GO TO 364                                                 191
  358 GO TO 370                                                              192
  360 KYL=1                                                                  193
  362 GO TO 366                                                              194
  364 KYL=104                                                                195
  366 TPC=RMARK                                                              196
  370 A(KYL) =TPC                                                            197
  386 IF(NCTR.GE.ILIM) GO TO 392                                             198
  388 NCTR=NCTR+1                                                            199
  390 GO TO 340                                                              200
C                                                                            201
  392 KOUT = 2                                                               202
  394 M= 10                                                                  203
C                                                                            204
  400 IF(XGL.AND.(MOD(LCTR,10).EQ.0))KQUIT=2                                 205
C                                                                            206
  600 WRITE (6,504)(A(I),I=1,104)                                            207
  602 IF(K864.NE.0) GO TO 620                                                208
C                                                                            209
  604 IF(MOD(LCTR,10))614,606,614                                            210
  606 XLABEL =FX+FLOAT(LCTR)*DX                                              211
  608 TEMP =ABS(XLABEL)                                                      212
  610 IF(TEMP.GE.1.E+7)LABOUT=2                                              213
  612 IF(TEMP.LT.1.E-7)XLABEL=0.                                             214
  613 WRITE(6,FLAB) XLABEL                                                   215
C                                                                            216
  614 LCTR = LCTR+1                                                          217
  616 GO TO (320,302),KQUIT                                                  218
C                                                                            219
  620 IF(.NOT.LS) GO TO 604                                                  220
  622 LS = .FALSE.                                                           221
  623 KT=K864-1                                                              222
  624 WRITE (6,FLS)(ELS(I),I=1,KT)                                           223
  626 GO TO 604                                                              224
  630 LABOUT= 4                                                              225
  700 GO TO (710,702,706,704,712),LABOUT                                     226
  702 WRITE (6,506)                                                          227
      GO TO 720                                                              228
C                                                                            229
  704 WRITE(6,520) NCTR                                                      230
  706 WRITE(6,508)(X(I),Y(I),I=1,2),K(1),K(2),K(3)                           231
  708 GO TO 720                                                              232
  712 WRITE(6,512)                                                           233
  710 WRITE(6,510)                                                           234
  720 RETURN                                                                 235
  500 FORMAT(2HPT)                                                           236
  502 FORMAT(2HPT,12X,3HX*E,I2)                                              237
  504 FORMAT(2HP ,26X,104A1)                                                 238
  506 FORMAT(2HPL,3X,10HBAD LABELS  )                                        239
  508 FORMAT(2HPL,5H N.G.,4G20.8,3I7)                                        240
  510 FORMAT(2HPL)                                                           241
  512 FORMAT(2HPL,16HERROR IN K ARRAY      )                                 242
  520 FORMAT(2HPL,18HX OUT OF ORDER. I=,I5)                                  243
      END                                                                    244
$IBFTC KHAR    ALTIO
       FUNCTION KHAR(XMAX)                                                     1
       KHAR = INT(ALOG(XMAX)/2.302585+40.0)-40                                 2
       RETURN                                                                  3
       END                                                                     4
$ORIGIN        ALPHA
$INCLUDE       INPUT
$ORIGIN        ALPHA
$INCLUDE       PRECAL,MHORIZ
$ORIGIN        ALPHA
$INCLUDE       COEF,COEFBB,HRB,AAK,BDRY12,BDRY34
$ORIGIN        ALPHA
$INCLUDE       SOR
$ORIGIN        ALPHA
$INCLUDE       PISTUG,PLOTMY,KHAR
$ORIGIN        BETA
$INCLUDE       SLAX,SLAV
$ORIGIN        ALPHA
$INCLUDE       TANG,SEARCH
$ORIGIN        BETA
$INCLUDE       VELOCY,VEL
      MIXED FLOW IMPELLER (NASA TN D-1186)
1.5       1000.     1000000.  1.        .003042   .0013516  796.0     0.
-84.88    -43.      .10555    -2.629    .05664    -.6649    .04891    -2.3434
   10   47   28   47   57   28    8   18
.000914   .001846   -80.0     -49.0     6.0
0.        .01214    .02651    .04766    .07360    0.
0.        -.6250    -1.2330   -1.8182   -2.2750   0.
.000914   .001846   -83.0     -41.5     6.0
0.        .00788    .02004    .04006    .06828    0.
0.        -.6310    -1.2310   -1.8206   -2.2954   0.
.001328   .001753   -60.5     -51.5     6.0
0.        .01307    .02552    .04172    .05280    0.
0.        -.1670    -.3370    -.5262    -.6269    0.
.001328   .001753   -63.0     -40.5     5.0
0.        .01073    .02493    .04172    0.
0.        -.2010    -.4070    -.5818    0.
-.03124   -.01514   .00025    .01065    .01853    .02651    .03460    .04281
.05115    .05964    .06828    .07709    .08607    .09524    .10461    .11417
.12726    .14073
.07586    .07662    .07874    .08091    .08294    .08531    .08802    .09108
.09447    .09820    .10228    .10674    .11146    .11656    .12200    .12778
.13602    .14487
.010533   .010045   .008724   .007420   .006316   .005354   .004532   .003831
.003235   .002728   .002299   .001936   .001629   .001370   .001151   .000979
.000825   .000724
    1    0    0    2    2    1    3
      MODIFIED TANDEM AXIAL TURBINE ROTOR
1.4       287.053   288.15    1.225     .03152    .011347             0.
48.0      -47.0     .02847    .021333   .02515    -.05459   .02441    -.003607
   10   32   29   49   58   20   76    2
.000762   .000381   50.0      -29.4     7.
          .00257    .00765    .01527    .02035    .02543
          .00925    .02118    .02988    .03020    .02643
.000762   .000381   25.0      -6.9      5.
          .00765    .02035    .02543
          .00714    .02039    .02094
.001778   .000381   -8.1      -48.8     4.0
0.        .00610    .01626    0.
0.        .00164    -.02463   0.
.001778   .000381   -19.7     -42.5     4.0
0.        .00610    .01372    0.
0.        -.01200   -.02745   0.
-1.       1.
.32385    .32385
.01       .01
    1    2    0    2    2    1    3
